diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 33961c0fb..96c4042ec 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,7 +20,7 @@ steps: - label: ":speedboat: GPU unit tests" commands: - - "$JULIA_PATH/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project -e 'using Pkg; Pkg.Registry.update(); Pkg.test()'" + - "$JULIA_PATH/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project --check-bounds=yes -e 'using Pkg; Pkg.Registry.update(); Pkg.test()'" depends_on: "init" notify: - github_commit_status: @@ -30,7 +30,7 @@ steps: env: CUDA_VISIBLE_DEVICES: "-1" commands: - - "$JULIA_PATH/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project -e 'using Pkg; Pkg.Registry.update(); Pkg.test()'" + - "$JULIA_PATH/julia-$JULIA_VERSION/bin/julia -O0 --color=yes --project --check-bounds=yes -e 'using Pkg; Pkg.Registry.update(); Pkg.test()'" depends_on: "init" notify: - github_commit_status: diff --git a/Manifest.toml b/Manifest.toml index f7baaf529..e18b111ac 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.2" manifest_format = "2.0" -project_hash = "4319b71589d68d891265adf31b9a9017f2b1b326" +project_hash = "37e11a3bdd396c973917441db34ded05b97faa85" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -362,6 +362,18 @@ git-tree-sha1 = "518ebd058c9895de468a8c255797b0c53fdb44dd" uuid = "61eb1bfa-7361-4325-ad38-22787b887f55" version = "0.26.5" +[[deps.GibbsSeaWater]] +deps = ["GibbsSeaWater_jll", "Libdl", "Test"] +git-tree-sha1 = "d1642ddc78d0754603d747d76fac0042a5a85104" +uuid = "9a22fb26-0b63-4589-b28e-8f9d0b5c3d05" +version = "0.1.3" + +[[deps.GibbsSeaWater_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c91ca76546871efaa1aefdd2b19cc41c3ead2160" +uuid = "6727f6b2-98ea-5d0a-8239-2f72283ddb11" +version = "3.5.2+0" + [[deps.Glob]] git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496" uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" diff --git a/Project.toml b/Project.toml index b1951ea73..3202b820c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,15 +1,17 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright and contributors"] -version = "0.10.5" +version = "0.11.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +GibbsSeaWater = "9a22fb26-0b63-4589-b28e-8f9d0b5c3d05" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" [compat] @@ -17,10 +19,12 @@ Adapt = "3, 4" CUDA = "5" Documenter = "1" EnsembleKalmanProcesses = "1" +GibbsSeaWater = "0.1" JLD2 = "0.4" KernelAbstractions = "0.9" Oceananigans = "0.91" Roots = "2" +SeawaterPolynomials = "0.3" StructArrays = "0.4, 0.5, 0.6" julia = "1.10" diff --git a/README.md b/README.md index 1dde63565..29712719e 100644 --- a/README.md +++ b/README.md @@ -130,8 +130,7 @@ grid = RectilinearGrid(GPU(), size = (256, 32), extent = (500meters, 100meters), ## Documentation -See the [documentation](https://oceanbiome.github.io/OceanBioME.jl) for full description of the software -package and more examples. +See the [documentation](https://oceanbiome.github.io/OceanBioME.jl) for full description of the software package and more examples, as well as full descriptions of the included models and parametrisations. ## Contributing If you're interested in contributing to the development of OceanBioME we would appreciate your help! diff --git a/docs/make.jl b/docs/make.jl index 66bf1f6dd..4a84eb77b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,10 +1,9 @@ using Documenter, DocumenterCitations, Literate using OceanBioME -using OceanBioME.SLatissimaModel: SLatissima -using OceanBioME.LOBSTERModel: LOBSTER -using OceanBioME.Boundaries.Sediments: SimpleMultiG, InstantRemineralisation -using OceanBioME.Boundaries: OCMIP_default, GasExchange +using OceanBioME: SLatissima, LOBSTER, NutrientPhytoplanktonZooplanktonDetritus +using OceanBioME.Sediments: SimpleMultiG, InstantRemineralisation +using OceanBioME: CarbonChemistry, GasExchange using Oceananigans.Grids: RectilinearGrid @@ -58,14 +57,14 @@ model_parameters = (LOBSTER(; grid = BoxModelGrid(), light_attenuation_model = n TwoBandPhotosyntheticallyActiveRadiation(; grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))), SimpleMultiG(; grid = BoxModelGrid()), InstantRemineralisation(; grid = BoxModelGrid()), - OCMIP_default, - GasExchange(; gas = :CO₂).condition.func, - GasExchange(; gas = :O₂).condition.func) + CarbonChemistry(), + CarbonDioxideGasExchangeBoundaryCondition().condition.func, + OxygenGasExchangeBoundaryCondition().condition.func) exchanged_gas(::Val{G}) where G = G model_name(model) = if Base.typename(typeof(model)).wrapper == GasExchange - "$(exchanged_gas(model.gas)) air-sea exchange" + ifelse(isa(model.water_concentration, CarbonChemistry), "CO₂", "O₂")*" air-sea exchange" else Base.typename(typeof(model)).wrapper end @@ -98,6 +97,7 @@ individuals_pages = [ component_pages = [ "Biogeochemical models" => bgc_pages, "Air-sea gas exchange" => "model_components/air-sea-gas.md", + "Carbon chemistry" => "model_components/carbon-chemistry.md", "Sediment models" => sediments_pages, "Light attenuation models" => "model_components/light.md", "Individuals" => individuals_pages, @@ -145,7 +145,7 @@ makedocs(sitename = "OceanBioME.jl", pages = pages, modules = [OceanBioME], plugins = [bib], - doctest = true, + doctest = false,#true, clean = true, checkdocs = :exports) diff --git a/docs/oceanbiome.bib b/docs/oceanbiome.bib index d08dc3643..e2c6412ee 100644 --- a/docs/oceanbiome.bib +++ b/docs/oceanbiome.bib @@ -293,3 +293,107 @@ @book{tan1998 publisher={Oxford University Press}, pages={271-291} } + +@article{dickson2007, + title = {Guide to best practices for ocean CO2 measurement}, + author = {Dickson, A.G. and Sabine, C.L and Christian, J.R.}, + year = {2007}, + doi = {10.25607/OBP-1342}, + journal = {North Pacific Marine Science Organization}, + pages = {191}, + volume = {PICES Special Publication 3; IOCCP Report 8} +} + +@article{branson2023, + author = {Branson, Oscar}, + doi = {10.5281/zenodo.7732476}, + month = mar, + journal = {Zenodo}, + title = {oscarbranson/cbsyst}, + url = {https://doi.org/10.5281/zenodo.7732476}, + year = 2023, +} + + +@article{humphreys2022, + author = {Humphreys, M. P. and Lewis, E. R. and Sharp, J. D. and Pierrot, D.}, + journal = {Geoscientific Model Development}, + number = {1}, + pages = {15--43}, + title = {PyCO2SYS v1.8: marine carbonate system calculations in Python}, + volume = {15}, + year = {2022} +} + + +@article{feistel2008, + abstract = {The specific Gibbs energy of seawater is determined from experimental data of heat capacities, freezing points, vapour pressures and mixing heats at atmospheric pressure in the range −6 to 80$\,^{\circ}$C in temperature and 0--120gkg--1 in absolute salinity. Combined with the pure-water properties available from the 1996 Release of the International Association for the Properties of Water and Steam (IAPWS-95), and the densities from the 2003 Gibbs function of seawater, a new saline part of the Gibbs function is developed for seawater that has an extended range of validity including elevated temperatures and salinities. In conjunction with the IAPWS 2006 Release on ice, the correct description of concentrated brines by the new formulations permits an accurate evaluation of sea ice properties down to salinity saturation temperatures. The new Gibbs function is expressed in terms of the temperature scale ITS-90. Its input variable for the concentration is absolute salinity, available from the new Reference-Composition Salinity Scale of 2008.}, + author = {Rainer Feistel}, + doi = {https://doi.org/10.1016/j.dsr.2008.07.004}, + issn = {0967-0637}, + journal = {Deep Sea Research Part I: Oceanographic Research Papers}, + number = {12}, + pages = {1639-1671}, + title = {A Gibbs function for seawater thermodynamics for −6 to 80$\,^{\circ}$C and salinity up to 120gkg--1}, + url = {https://www.sciencedirect.com/science/article/pii/S0967063708001489}, + volume = {55}, + year = {2008}, + bdsk-url-1 = {https://www.sciencedirect.com/science/article/pii/S0967063708001489}, + bdsk-url-2 = {https://doi.org/10.1016/j.dsr.2008.07.004} +} + + +@article{roquet2015, + author = {F. Roquet and G. Madec and Trevor J. McDougall and Paul M. Barker}, + journal = {Ocean Modelling}, + pages = {29-43}, + title = {Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard}, + volume = {90}, + year = {2015} +} + +@article{Ho2006, + abstract = {The SOLAS Air-Sea Gas Exchange (SAGE) Experiment was conducted in the western Pacific sector of the Southern Ocean. During SAGE, gas transfer velocities were determined using the 3He/SF6 dual gas tracer technique, and results were obtained at higher wind speeds (16.0 m s−1) than in previous open ocean dual tracer experiments. The results clearly reveal a quadratic relationship between wind speed and gas transfer velocity rather than a recently proposed cubic relationship. A new parameterization between wind speed and gas transfer velocity is proposed, which is consistent with previous 3He/SF6 dual tracer results from the coastal and open ocean obtained at lower wind speeds. This suggests that factors controlling air-sea gas exchange in this region are similar to those in other parts of the world ocean, and that the parameterization presented here should be applicable to the global ocean.}, + author = {Ho, David T. and Law, Cliff S. and Smith, Murray J. and Schlosser, Peter and Harvey, Mike and Hill, Peter}, + doi = {https://doi.org/10.1029/2006GL026817}, + eprint = {https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006GL026817}, + journal = {Geophysical Research Letters}, + number = {16}, + title = {Measurements of air-sea gas exchange at high wind speeds in the Southern Ocean: Implications for global parameterizations}, + url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006GL026817}, + volume = {33}, + year = {2006}, + bdsk-url-1 = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006GL026817}, + bdsk-url-2 = {https://doi.org/10.1029/2006GL026817} +} + +@article{Wanninkhof2014, + author = {Wanninkhof, Rik}, + title = {Relationship between wind speed and gas exchange over the ocean revisited}, + journal = {Limnology and Oceanography: Methods}, + volume = {12}, + number = {6}, + pages = {351-362}, + doi = {https://doi.org/10.4319/lom.2014.12.351}, + url = {https://aslopubs.onlinelibrary.wiley.com/doi/abs/10.4319/lom.2014.12.351}, + eprint = {https://aslopubs.onlinelibrary.wiley.com/doi/pdf/10.4319/lom.2014.12.351}, + year = {2014} +} + + +@article{Weiss1974, + author = {R.F. Weiss}, + doi = {https://doi.org/10.1016/0304-4203(74)90015-2}, + issn = {0304-4203}, + journal = {Marine Chemistry}, + number = {3}, + pages = {203-215}, + title = {Carbon dioxide in water and seawater: the solubility of a non-ideal gas}, + url = {https://www.sciencedirect.com/science/article/pii/0304420374900152}, + volume = {2}, + year = {1974}, + bdsk-url-1 = {https://www.sciencedirect.com/science/article/pii/0304420374900152}, + bdsk-url-2 = {https://doi.org/10.1016/0304-4203(74)90015-2} +} + + diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md index ec8fbe7c6..64b5ac576 100644 --- a/docs/src/appendix/library.md +++ b/docs/src/appendix/library.md @@ -5,7 +5,6 @@ Documenting the user interface. ## OceanBioME.jl ```@autodocs Modules = [OceanBioME] -private = false ``` ## Biogeochemical Models @@ -13,46 +12,47 @@ private = false ### Nutrient Phytoplankton Zooplankton Detritus ```@autodocs -Modules = [OceanBioME.NPZDModel] -private = false +Modules = [OceanBioME.Models.NPZDModel] ``` ### The Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) ```@autodocs -Modules = [OceanBioME.LOBSTERModel] -private = true +Modules = [OceanBioME.Models.LOBSTERModel] ``` ### Sugar kelp (Saccharina latissima) ```@autodocs -Modules = [OceanBioME.SLatissimaModel] -private = true +Modules = [OceanBioME.Models.SLatissimaModel] +``` + +### Carbon Chemistry + +```@autodocs +Modules = [OceanBioME.Models.CarbonChemistryModel] ``` ## Light Attenuation Models ```@autodocs Modules = [OceanBioME.Light] -private = false ``` -## Boundary Conditions +## Sediments ```@autodocs -Modules = [OceanBioME.Boundaries] -private = false +Modules = [OceanBioME.Models.Sediments] ``` +## Gas exchange boundary conditions + ```@autodocs -Modules = [OceanBioME.Boundaries.Sediments] -private = false +Modules = [OceanBioME.Models.GasExchangeModel, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity] ``` ## Box Model ```@autodocs Modules = [OceanBioME.BoxModels] -private = false ``` diff --git a/docs/src/model_components/air-sea-gas.md b/docs/src/model_components/air-sea-gas.md index 86c41900c..c8bdebd66 100644 --- a/docs/src/model_components/air-sea-gas.md +++ b/docs/src/model_components/air-sea-gas.md @@ -1,214 +1,61 @@ # [Air-sea gas exchange](@id air-sea-gas) -We currently have one air-sea gas exchange model implemented. The model, proposed by [Wanninkhof1992](@citet), calculates the solubility of the gas in the water dependent on the temperature and salinity, and calculates the flux depending on the solubility and mixing from the wind. - -Currently, the parameters for CO₂ and oxygen are included, but it would be very straightforward to add the parameters given in the original publication for other gases (e.g. inert tracers of other nutrients such as N₂). We also currently have a very simple formulation of the gas transfer velocity which depends on an average wind speed, but it would straightforwardly be expanded to permit variable wind speed (e.g. to simulate enhanced exchange from storms). - -It is straightforward to set up a boundary as an air-sea gas exchange: - -```@setup gasexchange -using OceanBioME -CO₂_flux = GasExchange(; gas = :CO₂) -using Oceananigans - -grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); - -model = NonhydrostaticModel(; grid, - biogeochemistry = LOBSTER(; grid, carbonates = true), - boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ), - tracers = (:T, :S)) +Air-sea gas transfer is typically parameterised as a function of temperature (``T``) and wind speed (``u_{10}``), and the concentration of the gas in the air (``C_a``) and in the surface water (``C_w``) in the form: +```math +F = k(u_{10}, T)(C_w - C_a), ``` +where `k` is the gas transfer velocity. +Our implementation is intended to be generic for any gas, so you can specify `air_concentration`, `water_concentration`, `transfer_velocity`, and `wind_speed` as any function in `GasExchange`, but we also provide constructors and default values for carbon dioxide and oxygen. +To setup carbon dioxide and/or oxygen boundary conditions you simply build the condition and then specify it in the model: ```@example gasexchange using OceanBioME - -CO₂_flux = GasExchange(; gas = :CO₂) -``` - -Where the symbol specifies the exchanged gas (currently `:CO₂` or `:O₂`). This can then be passed in the setup of a BGC model, for example: - -```@example gasexchange +CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition() +O₂_flux = OxygenGasExchangeBoundaryCondition() using Oceananigans grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); model = NonhydrostaticModel(; grid, - biogeochemistry = LOBSTER(; grid, carbonates = true), - boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ), + biogeochemistry = LOBSTER(; grid, carbonates = true, oxygen = true), + boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), + O₂ = FieldBoundaryConditions(top = O₂_flux)), tracers = (:T, :S)) ``` -## Model equations - -### Gas transfer - -The gas flux is given by: - -```math -F = k(C_w - \alpha C_a), -``` - -where ``C_w`` is the concentration in the water, ``C_a`` the concentration in the air, ``\alpha`` the Oswald solubility coefficient, and ``k`` the gas transfer velocity. For carbon dioxide the flux is modified to: - -```math -F = k\beta\rho_o(pCO_{2w} - pCO_{2a}), -``` - -where ``pCO_{2w}`` and ``pCO_{2a}`` are the partial pressure of carbon dioxide in the water and air, ``\beta`` is the Bunsen Solubility Coefficient, and ``\rho_o`` is the density of the water. +!!! note -The gas transfer velocity is parameterised by the wind speed and Schmidt number, which in turn is parameterised by the temperature and salinity. The gas transfer velocity is given by: + All gas exchange models require temperature (`T`) to be present in the model, and carbon dioxide requires sailinity (`S`), total inorganic carbon (`DIC`), and alkalinity (`Alk`), and optionally can take silicate and phosphate where there names are specified in the keyword argument `silicate_and_phosphate_names` -```math -k = 1.08\times10^{-6}u^2\left(\frac{Sc}{660}\right)^{-1/2}, -``` - -where ``u`` is the winds speed 10m above the surface, and Sc is the Schmidt number parameterised as: - -```math -Sc = A - BT + CT^2 - DT^3, -``` - -where ``T`` is temperature in Kelvin and the other parameters are dependent on the gas type and given in [Parameters](@ref parameters). - -The solubilities are given by: - -```math -\alpha = 0.00367 T \exp{A_1 + 100\frac{A_2}{T} + A_3 \ln{\frac{T}{100}} + S\left(B_1 + \frac{B_2}{T} + \frac{B_3}{T^2}\right)}, -``` - -and - -```math -\beta = \exp{A_1 + 100\frac{A_2}{T} + A_3 \ln{\frac{T}{100}} + S\left(B_1 + \frac{B_2}{T} + \frac{B_3}{T^2}\right)}, -``` - -where ``S`` is salinity in practical units and the other default parameters are given in [Parameters](@ref parameters). - -### Partial pressure of carbon dioxide - -We currently do not have the full OCMIP partial pressure formulation, instead we follow the simplified formulation (as used in [Aumont2015](@citet)) where the partial pressure of CO``_2`` (``\mu``atm) in gas is found from: - -```math -pCO_{2a} = f_{CO_2}P_a, -``` - -where ``f_{CO_2}`` is the air fraction (ppmv), and ``P_a`` is the atmospheric pressure (atm). The ocean partial pressure (``\mu``atm) is derived from the dissolved inorganic carbon content, and alkalinity, as described by [tan1998](@citet), from the equilibrium of the following chemical system: - -```math -\ce{CO_2(g)} \ce{<=> CO_2(aq)},\ K_0=\frac{[\ce{O_2(aq)}]}{\ce{pCO_2}}, -``` - -```math -\ce{CO_2(aq) + H_2O}\ce{<=> H^+ + HCO^-_3},\ K_1=\frac{\ce{[HCO^-_3][H^+]}}{\ce{[CO_2(aq)]}}, -``` - -```math -\ce{HCO^-_3} \ce{<=> H^+ + CO^{2-}_3},\ K_2=\frac{\ce{[CO^{2-}_3][H^+]}}{\ce{HCO^-_3}}, -``` - -```math -\ce{H_2O} \ce{<=> H^+ + OH^-}. -``` - -These have rates constants: - -```math -K_0 = \frac{[\ce{CO_2(aq)}]}{\ce{pCO_2}}, -``` - -```math -K_1 = \frac{\ce{[HCO^-_3][H^+]}}{\ce{[CO_2(aq)]}}, -``` - -```math -K_2=\frac{\ce{[CO^{2-}_3][H^+]}}{\ce{HCO^-_3}}, -``` - -```math -K_w = \ce{[OH^-][H^+]}. -``` - -In this system DIC and Alk are defined to be: - -```math -\ce{DIC} = \ce{[CO_2(aq)] + [HCO^-_3] + [CO_3^{2-}]}, -``` - -```math -\ce{Alk} = \ce{[HCO^-_3] + 2[CO^{2-}_3] + [B(OH)^-_4] + [OH^-] - [H^+] + [OH^-] \pm [minor species]}. -``` - -To solve this we must find the Boric acid dissociation from: - -```math -\ce{B(OH)_3 <=>H^+ + B(OH)^-_4},\ K_B = \frac{\ce{[B(OH)^-_4][H^+]}}{\ce{[B(OH)_3]}}, -``` - -resulting in the equilibrium constant: - -```math -K_B = \frac{\ce{[B(OH)^-_4][H^+]}}{\ce{[B(OH)_3]}}. -``` - -Finally, taking the DIC and Alkalinity in micro equivalents (i.e. scaled by ``10^{-3}/\rho_o`` from mmol C/m``^3``) denoted by ``\bar{DIC}`` and ``\bar{Alk}``, the carbonate alkalinity is given by ``Alk_C = \bar{Alk} + [B(OH)^-_4] + [OH^-]``, and we define the boron content, ``B``, to be ``R_{B:S}S`` mol/kg. +## Model equations -The system can be rearranged to give two equations for the carbonate alkalinity: +### Gas transfer velocity +The default gas transfer velocity (`ScaledTransferVelocity`) returns a velocity in the form: ```math -\ce{Alk_C} = DIC\frac{K_1H + 2K^1K^2}{H^2 + K_1H+K_1K_2}, +k(u_{10}, T) = cu_{10}^2\left(\frac{Sc(T)}{660}\right)^{-1/2}, ``` - -```math -\ce{Alk_C} = \ce{Alk} - \frac{K_W}{H} - \frac{K_BB}{K_B + H}. -``` - -The equilibrium constants (``K_0``, ``K_1``, and ``K_2``) are given by: +where ``c`` is a coefficient (`coeff`) which typically is wind product specific with default value ``0.266`` cm/hour from [Ho2006](@citet), and ``Sc`` is gas specific the temperature dependent Schmidt number (the dimensionless ratio of momentum and mass diffusivity) specified as `schmidt_number` which can be any function of temperature. The default parameterisations is the 4th order polynomial formulation of [Wanninkhof2014](@citet). -```math -K_0 = \exp{\left(A_0 + \frac{B_0}{T} + C_0\log(\frac{T}{100}) + S (D_0 - E_0T + F_0\left(\frac{T}{100}\right)^2)\right)}, -``` +Currently, the parameters for CO₂ and oxygen are included, but it would be very straightforward to add the parameters given in the original publication for other gases (e.g. inert tracers of other nutrients such as N₂). -```math -K_1 = \exp{\left(A_1 - \frac{B_1}{T} - C_1 \log(T) - (D_1 + \frac{E_1}{T})\sqrt{S} + F_1S - G_1S^{1.5}\right)}, -``` +### Carbon dioxide partial pressure -```math -K_2 = \exp{\left(A_2 - \frac{B_2}{T} - C_2 \log(T) - (D_2 + \frac{E_2}{T})\sqrt{S} + F_2S - G_1S^{1.5}\right)}, -``` +For most gasses the water concentration `C_w` is simply taken directly from the biogeochemical model or another tracer (in which case `water_concentration` should be set to `TracerConcentration(:tracer_name)`), but for carbon dioxide the fugacity (``fCO_2``) must be derived from the dissolved inorganic carbon (`DIC`) and `Alk`alinity by a `CarbonChemistry` model (please see the docs for [CarbonChemistry](@ref carbon-chemistry)), and used to calculate the partial pressure (``pCO_2``). +The default parameterisation for the partial pressure (`CarbonDioxideConcentration`) is given by [dickson2007](@citet) and defines the partial pressure to be the mole fraction ``x(CO_2)`` multiplied by the pressure, ``P``, related to the fugacity by: ```math -K_b = \exp{\left(\frac{A_b - B_b\sqrt{S} - C_bS + D_b S^{1.5} - E_bS^2}{T} + F_b + G_B \sqrt{S} + H_b S - (I_b + J_b \sqrt{S} + S) \log(T) + L_b \sqrt{S} T\right)}, +fCO_2 = x(CO_2)P\exp\left(\frac{1}{RT}\int_0^P\left(V(CO_2)-\frac{RT}{P'}\right)dP'\right). ``` - +The volume (``V``) is related to the gas pressure by the virial expression: ```math -K_w = \exp{\left(A_w + \frac{B_w}{T} + C_w \ln T + \sqrt{S}\left(D_w + \frac{E_w}{T} + F_w \ln T\right)G_w S\right)}. +\frac{PV(CO_2)}{RT}\approx1+\frac{B(x, T)}{V(CO_2)}+\mathcal{O}(V(CO_2)^{-2}), ``` - -We solve these equations iteratively from an initial guess of ``pH=8`` to find ``H``, from which the partial pressure of ``CO_2`` is calculated as: - +and the first virial coefficient ``B`` for carbon dioxide in air can be approximated as: ```math -pCO_2 = 10^6\frac{Alk_CH^2}{K_0(K_1H + 2 K_1K_2)} +B_{CO_2-\text{air}} \approx B_{CO_2}(T) + 2x(CO_2)\delta_{CO_2-\text{air}}(T), ``` +where ``\delta`` is the cross virial coefficient. -### Parameter variable names - -Gas transfer - -| Symbol | Variable name | Units | -|------------------------------------------|----------------------------|-------------| -| ``A``, ``B``, ``C``, ``D`` | `schmidt_params` | - | -| ``A_i`` and ``B_i`` for ``i\in[1, 2, 3]``| `solubility_params` | - | - -Partial pressure of carbon dioxide - -| Symbol | Variable name | Units | -|------------------------------------------|----------------------------|-------------| -| ``L_0`` where ``L`` is a Latin character | `solubility` | - | -| ``L_1`` | `bicarbonate_dissociation` | - | -| ``L_2`` | `carbonate_dissociation` | - | -| ``L_b`` | `boric_acid_dissociation` | - | -| ``L_w`` | `water_dissociaiton` | - | -| ``R_{B:S}`` | `boron_ratio` | mol B / psu | -| ``\alpha`` | `thermal_expansion` | 1 / °K | -| ``\beta`` | `haline_contraction` | 1 / psu | \ No newline at end of file +``B_{CO_2}`` and ``\delta_{CO_2-\text{air}}`` are parameterised by [Weiss1974](@citet) and reccomended in [dickson2007](@citet) as fourth and first order polynomials respectively. diff --git a/docs/src/model_components/carbon-chemistry.md b/docs/src/model_components/carbon-chemistry.md new file mode 100644 index 000000000..e9d94c61d --- /dev/null +++ b/docs/src/model_components/carbon-chemistry.md @@ -0,0 +1,354 @@ +# [Carbon Chemistry](@id carbon-chemistry) + +Our carbon chemistry model follows the best practice described by [dickson2007](@citet), as implemented by e.g. `cbsyst` [branson2023](@citep) and `CO2SYS` [humphreys2022](@citep) in Python. + +The carbon chemistry model is primarily used for diagnosing the partial pressure of carbon dioxide in the surface water for [gas exchange with the air](@ref air-sea-gas), but is capable of diagnosing other species such as carbonate concentration, useful in the calculation of calcite dissolution. +The model works by computing the pH from the total dissolved inorganic carbon and total alkalinity (and optionally silicate, phosphate, boron, sulfate, and fluoride concentration), along with temperature and salinity, which can then be used to find the concentration of different species. +We will first describe how to use the model, followed by the underlying chemistry and parameterisations. + +## Using the model +To use the carbon chemistry model it is constructed by simply writing: +```@example carbon-chem +using OceanBioME + +carbon_chemistry = CarbonChemistry() +``` + +### Computing ``fCO_2`` and ``pH`` +To compute the fugacity of carbon dioxide (``fCO_2``) you call the model with the `DIC`, `Alk`alinity, `T`emperature, and `S`alinity: +```@example carbon-chem +DIC = 2145 +Alk = 2448 +T = 25.4 +S = 36.45 +carbon_chemistry(; DIC, Alk, T, S) +``` +This is sufficient when computing ``fCO_2`` at the surface, but if we wanted to know ``fCO_2`` at depth where there is higher pressure we can specify the pressure in bar like: +```@example carbon-chem +carbon_chemistry(; DIC, Alk, T, S, P = 100) +``` + +We may also be interested in the pH so we can request that be outputted: +```@example carbon-chem +carbon_chemistry(; DIC, Alk, T, S, return_pH = true) +``` + +These function calls assume a constant boron, sulfate, and fluoride ratio relative to the salinity (as described below), but can be specified instead: +```@example carbon-chem +carbon_chemistry(; DIC, Alk, T, S, boron = 0) +``` + +And the silicate and phosphate concentrations are assumed to be zero but can similarly be specified: +```@example carbon-chem +carbon_chemistry(; DIC, Alk, T, S, silicate = 2, phosphate = 1) +``` + +The same code can also be used to compute ``fCO_2`` when the pH is already known by passing it in the same way: +```@example carbon-chem +carbon_chemistry(; DIC, pH = 8.1, T, S) +``` + +Finally, for slightly improved accuracy you may wish to specify that the seawater density calculation, is an intermediary step in the calculations above, is computed using the full TEOS-10 [feistel2008](@citet) standard by setting this as the `density_function` when you construct the carbon chemistry model: +```@example carbon-chem +using OceanBioME.Models: teos10_density +carbon_chemistry = CarbonChemistry(; density_function = teos10_density) +``` +But this comes at the cost of significantly increased computational expense and GPU incompatability. + +During the conversion from practical to absolute salinity units the location can then also be entered for (very slightly, in this example ~``10^{-4}\mu``atm) improved accuracy: +```@example carbon-chem +carbon_chemistry(; DIC, Alk, T, S, lon = -31.52, lat = 33.75) +``` + +The default uses the polynomial approximation described in [roquet2015](@citet) as provided by [`SeawaterPolynomials.jl`](https://github.com/CliMA/SeawaterPolynomials.jl/). + +### Computing the carbonate concentration +So that this model can be used in calcite dissolution models it can also return the carbonate saturation by calling the function `carbonate_saturation` +```@example carbon-chem +using OceanBioME.Models.CarbonChemistryModel: carbonate_saturation + +carbonate_saturation(carbon_chemistry; DIC, Alk, T, S) +``` +This function takes all of the same arguments (e.g. `boron`) as `carbon_chemistry` above. + +## Chemistry +### pH computation +When carbon dioxide is dissolved in seawater it dissociates into carbonate and bicarbonate species according to the equilibria: + +```math +\ce{CO_2(g)} \ce{<=> CO_2(aq)}, +``` + +```math +\ce{CO_2(aq)} + \ce{H_2O} \ce{<=> H_2CO_3(aq)}, +``` + +```math +\ce{CO_2(aq) + H_2O}\ce{<=> H^+ HCO^-_3}, +``` + +```math +\ce{HCO^-_3} \ce{<=> H^+ + CO^{2-}_3}, +``` + +from which we define the total dissolved inorganic carbon (`DIC`) to be ``DIC = [\ce{CO_2(aq)}] + [\ce{HCO^-_3}] + [\ce{CO^{2-}_3}]``. +The equilibrium constants for these equations are defined as: + +```math +K_0=\frac{[\ce{CO_2(aq)}]}{\ce{pCO_2}}, +``` + +```math +K_1=\frac{\ce{[HCO^-_3][H^+]}}{\ce{[CO_2(aq)]}}, +``` +and +```math +K_2=\frac{\ce{[CO^{2-}_3][H^+]}}{\ce{HCO^-_3}}. +``` + +These equilibria depend on the total hydrogen ion concentration ``[H^+]``, which depends on the concentration of acids and buffering of bases. + +#### Alkalinity: acid-base buffering + +Alkalinity is a measure of the buffering capacity of seawater and we can approximate it as: +```math +Alk = [\ce{HCO_3^-}] + 2[\ce{CO_3^{2-}}] + [\ce{B(OH)_4^-}] + [\ce{OH^-}] + [\ce{HPO_4^{2-}}] + 2[\ce{PO_4^{3-}}] + [\ce{SiO(OH)_3^-}] + [\ce{NH_3}] + [\ce{HS^-}] - [\ce{H^+}] - [\ce{HSO_4^-}] - [\ce{HF}] - [\ce{H_3PO_4}] + \text{ minor acids and bases}, +``` +where the "minor" species are either in sufficiently low concentrations as to be neglected (sometimes the photphate, silicate and sulfate species can also be neglected), and we shall neglect ``[\ce{NH_3}]``. + +Each of these acid and base species are also in an equilibrated dissociation state given below. + +#### Hydrogen sulfate +Sulfate in the form of hydrogen sulfate dissociates to form sulfate ions in the equilibrium +```math +\ce{HSO_4^-}\ce{<=> H^+ + SO_4^{2-}}, +``` +with an equilibrium constant given by +```math +K_S=\frac{\ce{[SO_4^{2-}][H^+]}}{\ce{HSO_4^-}}. +``` + +#### Boric acid +Boron in the form of boric acid equilibrates with water in the reaction +```math +\ce{B(OH)_3+H_2O}\ce{<=> H^+ + B(OH)_4^-}, +``` +with equilibrium constant +```math +K_B=\frac{\ce{[B(OH)_4^-][H^+]}}{\ce{B(OH)_3}}. +``` + +#### Hydrogen fluoride +Hydrogen fluoride dissociated in the equilibrium +```math +\ce{HF}\ce{<=> H^+ + F^-}, +``` +with equilibrium constant +```math +K_F=\frac{\ce{[F^-][H^+]}}{\ce{HF}}. +``` + +#### Phosphoric acid +Phosphoric acid undergoes a three stage dissociation given by the equilibrium +```math +\ce{H_3PO_4}\ce{<=> H^+ + H_2PO_4^-}, +``` + +```math +\ce{H_2PO_4^-}\ce{<=> H^+ + HPO_4^{2-}}, +``` + +```math +\ce{HPO_4^{2-}}\ce{<=> H^+ + PO_4^{3-}}, +``` + +with equilibrium constants + +```math +K_{P1} = \frac{\ce{[H^+][H_2PO_4^-]}}{\ce{H_3PO_4}}, +``` + +```math +K_{P2} = \frac{\ce{[H^+][HPO_4^{2-}]}}{\ce{H_2PO_4^-}}, +``` + +```math +K_{P3} = \frac{\ce{[H^+][PO_4^{3-}]}}{\ce{HPO_4^{2-}}}. +``` + +#### Silicic acid +Silicic acid dissociates in the equilibrium +```math +\ce{Si(OH)_4}\ce{<=> H^+ + SiO(OH)_3^-}, +``` + +with equilibrium constant +```math +K_{Si} = \frac{\ce{[H^+][SiO(OH)_3^-]}}{\ce{Si(OH)_4}}. +``` + +#### Water +Finally, water dissociates in the equilibrium +```math +\ce{H_2O}\ce{<=> H^+ + OH^-}, +``` + +with equilibrium constant +```math +K_w = \\ce{[H^+][OH^-]}. +``` + +#### Alkalinity equilibration +From these rate constants we can rewrite the total alkalinity (given above) in terms of only the rate constants, total hydrogen ion concentration (``\ce{H^+}``), the total dissolved inorganic carbon (``[DIC]``), boron (``[\text{B}]``), phosphate (``[\text{P}]``), silicate (``[\text{Si}]``), Sulfate (``[\text{Sulf}]``), and fluorine (``[\text{F}]``) content of the water, by rearranging the equations above. +This results in a form of the total alkalinity: +```math +\begin{align} +Alk &\approx \frac{[DIC][\ce{H^+}]K_1}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2}\\ + &+ \frac{2[DIC]K_1K_2}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2}\\ + &+ \frac{[\text{B}]}{1+[\ce{H^+}]/K_B}\\ + &+ \frac{K_w}{[\ce{H^+}]}\\ + &+ \frac{[\text{P}][\ce{H^+}]K_{P1}K_{P2}}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}\\ + &+ \frac{2[\text{P}]K_{P1}K_{P2}K_{P3}}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}\\ + &+ \frac{[\text{Si}]}{1+[\ce{H^+}]/K_{Si}}\\ + &- \frac{[\ce{H^+}]}{1+\text{S}/K_S}\\ + &- \frac{[\text{Sulf}]}{1+K_S/[\ce{H^+}]/(1+S/K_S)}\\ + &- \frac{[\text{F}]}{1+K_F/[\ce{H^+}]}\\ + &- \frac{[\ce{H^+}]^2}{[\ce{H^+}]^3+K_{P1}[\ce{H^+}]^2+K_{P1}K_{P2}[\ce{H^+}] + K_{P1}K_{P2}K_{P3}}. +\end{align} +``` + +This gives us a second equation from total alkalinity (which we must already know) with one unknown, ``[\ce{H^+}]``. +Our model solves for ``[\ce{H^+}]`` using a bisection method to an accuracy of ``10^{-10}``, having approximated the equilibrium constants from parametrisations described below. +We can then determine the pH as, ``pH = -\log_{10}([\ce{H^+}])``. + +### Carbon dioxide +Now that we know ``[DIC]``, ``Alk``, and ``[\ce{H^+}]`` we can return to the equation for total dissolved inorganic carbon to find the concentration of aqueous carbon dioxide +```math +[CO_2(aq)] = \frac{[DIC][\ce{H^+}]^2}{[\ce{H^+}]^2 + K_1[\ce{H^+}] + K_1K_2}, +``` + +from which we determine gas phase carbon dioxide concentration as +```math +fCO_2 = \frac{[CO_2(aq)]}{K_0}, +``` +in atmospheres. + +### Carbonate concentration and calcite saturation +Similarly we can also diagnose the calcite concentration +```math +[CO_3^{2-}] = \frac{[DIC]K_1K_2}{[\ce{H^+}]([\ce{H^+}] + K_1)+K_1K_2}. +``` + +This concentration is important in the dissolution of calcium carbonate which reacts according to the equilibrium +```math +\ce{CaCO_3(aq)}\ce{<=> Ca^+ + CO_3^{2-}}, +``` +which has an equilibrium constant +```math +K_{SP} = [\ce{Ca^{2+}}][\ce{CO_3^{2-}}]. +``` + +The calcite saturation can then be defined as ``\Omega=\frac{[\ce{CO_3^{2-}}]}{[\ce{CO_{3, saturation}^{2-}}]}`` which can be diagnosed as: +```math +\Omega = \frac{[\ce{Ca^{2+}}][\ce{CO_3^{2-}}]}{K_{SP}}. +``` + +### Missing pieces +In most cases the chemistry described above requires information about more elements that is usually available. +This means that we must parameterise their concentrations, usually this results in the boron, sulfate, and fluoride concentrations being set as constant ratios of the salinity. +Usually these ratios are: +```math +\begin{align} +[\text{B}] &= \frac{0.000232}{10.811}\frac{S}{1.80655},\\ +[\text{S}] &= \frac{0.14}{96.06} \frac{S}{1.80655},\\ +[\text{F}] &= \frac{0.000067}{18.9984} \frac{S}{1.80655}.\\ +\end{align} +``` +We use these ratios by default in the model (but they can be changed when calling the models functions described below). +Additionally, silicate and phosphate concentrations are often unavailable (both in observations and models), so are by default set to zero and their alkalinity contribution assumed to be small. + +The "ionic strength" must also be parameterised for some equilibrium constants and is usually assumed to be: +```math +Is = \frac{19.924}{1000 - 1.005S}, +``` +but these parameters are also variable in the model. + +## Model parameterisation +The chemical system described above has a large number of equilibrium constants, the constants are typically assumed to depend on the temperature, salinity, and pressure (hydrostatic pressure from depth underwater). +[//]: # (This sentence is not currently actually true, but will be by the time of PR merge) +By default, this model parameterises them based on the "best practice" guidelines of [dickson2007](@citet). +These parameterisations are: +```@example carbon-chem +using OceanBioME.Models.CarbonChemistryModel: K0, K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi, KSP_calcite, KSP_aragonite +K0() # Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values) +``` +```@example carbon-chem +K1() # Millero (1995, Geochim. Cosmochim. Acta, 59, 664) +``` +```@example carbon-chem +K2() # Millero (1995, Geochim. Cosmochim. Acta, 59, 664) +``` +```@example carbon-chem +KB() # Dickson (1990, Deep-Sea Res., 37, 755–766) +``` +```@example carbon-chem +KW() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) +``` +```@example carbon-chem +KS() # Dickson (1990, Chem. Thermodyn., 22, 113–127) +``` +```@example carbon-chem +KF() # Dickson and Riley (1979, Mar. Chem., 7, 89–99) +``` +```@example carbon-chem +KP1() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) +``` +```@example carbon-chem +KP2() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) +``` +```@example carbon-chem +KP3() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) +``` +```@example carbon-chem +KSi() # Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) +``` +```@example carbon-chem +KSP_calcite() # Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341) +``` +```@example carbon-chem +KSP_aragonite() # Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341) +``` +The pressure corrections from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341) are: +```@example carbon-chem +K1().pressure_correction +``` +```@example carbon-chem +K2().pressure_correction +``` +```@example carbon-chem +KB().pressure_correction +``` +```@example carbon-chem +KW().pressure_correction +``` +```@example carbon-chem +KS().pressure_correction +``` +```@example carbon-chem +KF().pressure_correction +``` +```@example carbon-chem +KP1().pressure_correction +``` +```@example carbon-chem +KP2().pressure_correction +``` +```@example carbon-chem +KP3().pressure_correction +``` +```@example carbon-chem +KSP_calcite().pressure_correction +``` +```@example carbon-chem +KSP_aragonite().pressure_correction +``` \ No newline at end of file diff --git a/docs/src/model_components/sediments/simple_multi_g.md b/docs/src/model_components/sediments/simple_multi_g.md index f910ab286..5c1a26696 100644 --- a/docs/src/model_components/sediments/simple_multi_g.md +++ b/docs/src/model_components/sediments/simple_multi_g.md @@ -4,7 +4,7 @@ This model, proposed by [Soetaert2000](@citet), is a "G class" model that evolve It is straightforward to set up: -```jldoctest simplemultig; filter = r".*@ OceanBioME.Boundaries.Sediments.*" +```jldoctest simplemultig; filter = r".*@ OceanBioME.Models.Sediments.*" using OceanBioME, Oceananigans, OceanBioME.Sediments grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)) @@ -13,7 +13,7 @@ sediment_model = SimpleMultiG(; grid) # output ┌ Warning: Sediment models are an experimental feature and have not yet been validated. -└ @ OceanBioME.Boundaries.Sediments ~/OceanBioME.jl/src/Boundaries/Sediments/simple_multi_G.jl:102 +└ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:104 [ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC. Single-layer multi-G sediment model (Float64) ``` diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 231039724..ccb9290e3 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -113,35 +113,31 @@ Hence, we define the nutrient forcing using as the negative of the phytoplankton @inline (bgc::NutrientPhytoplankton)(::Val{:N}, args...) = -bgc(Val(:P), args...) ``` -Finally, we need to define some functions to allow us to update the time-dependent parameters (in this case the PAR and temperature, ``T``): - -```@example implementing -using OceanBioME: BoxModel -import OceanBioME.BoxModels: update_boxmodel_state! - -function update_boxmodel_state!(model::BoxModel{<:NutrientPhytoplankton, <:Any, <:Any, <:Any, <:Any, <:Any}) - getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time) - getproperty(model.values, :T) .= model.forcing.T(model.clock.time) -end -``` - Now we can run an example similar to the [LOBSTER box model example](@ref box_example): ```@example implementing using OceanBioME, Oceananigans.Units +using Oceananigans.Fields: FunctionField const year = years = 365days @inline PAR⁰(t) = 500 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2 +clock = Clock(; time = 0.0) + z = -10 # specify the nominal depth of the box for the PAR profile -@inline PAR(t) = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay +@inline PAR_func(t) = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay + +PAR = FunctionField{Center, Center, Center}(PAR_func, BoxModelGrid(); clock) @inline temp(t) = 2.4 * cos(t * 2π / year + 50days) + 26 -biogeochemistry = NutrientPhytoplankton() +biogeochemistry = Biogeochemistry(NutrientPhytoplankton(); + light_attenuation = PrescribedPhotosyntheticallyActiveRadiation(PAR)) -model = BoxModel(; biogeochemistry, forcing = (; PAR, T = temp)) +model = BoxModel(; biogeochemistry, + prescribed_tracers = (; T = temp), + clock) set!(model, N = 15, P = 15) @@ -174,7 +170,7 @@ axP = Axis(fig[1, 2], ylabel = "Phytoplankton \n(mmol N / m³)") lines!(axP, times / year, P[1, 1, 1, :], linewidth = 3) axPAR= Axis(fig[2, 1], ylabel = "PAR (einstein / m² / s)", xlabel = "Time (years)") -lines!(axPAR, times / year, PAR.(times), linewidth = 3) +lines!(axPAR, times / year, PAR_func.(times), linewidth = 3) axT = Axis(fig[2, 2], ylabel = "Temperature (°C)", xlabel = "Time (years)") lines!(axT, times / year, temp.(times), linewidth = 3) @@ -200,9 +196,9 @@ biogeochemical_drift_velocity(bgc::NutrientPhytoplankton, ::Val{:P}) = Another aspect that OceanBioME includes is sediment models. Doing this varies between sediment models, but for the most generic and simplest, all we need to do is add methods to two functions: ```@example implementing -using OceanBioME.Boundaries.Sediments: sinking_flux +using OceanBioME.Sediments: sinking_flux -import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers +import OceanBioME.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers @inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) @@ -231,7 +227,7 @@ using OceanBioME.Sediments: InstantRemineralisation # define the grid -grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = (32, ), extent = (100, )) +grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = (32, ), x = 1, y = 1, z = (-100, 0)) # setup the biogeochemical model diff --git a/examples/column.jl b/examples/column.jl index 2ddaff8d9..3f335a0f8 100644 --- a/examples/column.jl +++ b/examples/column.jl @@ -14,7 +14,6 @@ # ## Model setup # We load the packages and choose the default LOBSTER parameter set using OceanBioME, Oceananigans, Printf -using OceanBioME.SLatissimaModel: SLatissima using Oceananigans.Fields: FunctionField, ConstantField using Oceananigans.Units @@ -52,7 +51,7 @@ biogeochemistry = LOBSTER(; grid, carbonates = true, scale_negatives = true) -CO₂_flux = GasExchange(; gas = :CO₂) +CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition() clock = Clock(; time = 0.0) T = FunctionField{Center, Center, Center}(temp, grid; clock) @@ -116,10 +115,12 @@ carbon_export = zeros(length(times)) using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity -for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) +for (n, t) in enumerate(times) + clock.time = t + + air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = DIC[n], Alk = Alk[n], T, S)) + carbon_export[n] = (sPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/data_forced.jl b/examples/data_forced.jl index 409580655..4e748b655 100644 --- a/examples/data_forced.jl +++ b/examples/data_forced.jl @@ -84,7 +84,7 @@ biogeochemistry = LOBSTER(; grid, carbonates = true, scale_negatives = true) -CO₂_flux = GasExchange(; gas = :CO₂) +CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition() T = FunctionField{Center, Center, Center}(t_function, grid; clock) S = FunctionField{Center, Center, Center}(s_function, grid; clock) @@ -155,10 +155,13 @@ carbon_export = zeros(length(times)) using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity -for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], t_function(1, 1, 0, t), s_function(1, 1, 0, t)) - carbon_export[i] = (sPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + - bPOM[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) +for (n, t) in enumerate(times) + clock.time = t + + air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = DIC[n], Alk = Alk[n], T, S)) + + carbon_export[n] = (sPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:sPOM)).w[1, 1, grid.Nz-20] + + bPOM[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:bPOM)).w[1, 1, grid.Nz-20]) * redfield(Val(:sPOM), model.biogeochemistry) end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/examples/eady.jl b/examples/eady.jl index 8395a9c53..67cafc98b 100644 --- a/examples/eady.jl +++ b/examples/eady.jl @@ -57,7 +57,7 @@ biogeochemistry = LOBSTER(; grid, carbonates = true, open_bottom = true) -DIC_bcs = FieldBoundaryConditions(top = GasExchange(; gas = :CO₂)) +DIC_bcs = FieldBoundaryConditions(top = CarbonDioxideGasExchangeBoundaryCondition()) # Model instantiation model = NonhydrostaticModel(; grid, diff --git a/examples/kelp.jl b/examples/kelp.jl index 8157e082e..2d1921824 100644 --- a/examples/kelp.jl +++ b/examples/kelp.jl @@ -48,7 +48,7 @@ Lx, Ly = 20meters, 20meters grid = RectilinearGrid(architecture, size=(1, 1, 50), extent=(Lx, Ly, 200)) # Specify the boundary conditions for DIC and O₂ based on the air-sea CO₂ and O₂ flux -CO₂_flux = GasExchange(; gas = :CO₂) +CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition() clock = Clock(; time = 0.0) T = FunctionField{Center, Center, Center}(temp, grid; clock) @@ -139,10 +139,13 @@ carbon_export = zeros(length(times)) using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity -for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) - carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + - bPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] +for (n, t) in enumerate(times) + clock.time = t + + air_sea_CO₂_flux[n] = CO₂_flux.condition.func(1, 1, grid, clock, (; DIC = DIC[n], Alk = Alk[n], T, S)) + + carbon_export[n] = sPOC[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + + bPOC[n][1, 1, grid.Nz-20] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] end # Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. diff --git a/src/Boundaries/Boundaries.jl b/src/Boundaries/Boundaries.jl deleted file mode 100644 index d19f17aef..000000000 --- a/src/Boundaries/Boundaries.jl +++ /dev/null @@ -1,40 +0,0 @@ -""" -Boundary conditions for air/sea and sediment flux. - -Currently implemented: - -- gasexchange [Wanninkhof1992](@cite) - - Generic air-sea flux model described by [Wanninkhof1992](@citet) but only setup for CO₂ and O₂. - - Forces the DIC and oxygen fields, and requires temp (in centigrade) and salinity, plus - current DIC and ALK concentration. - -- Sediments - - Soetaert [Soetaert2000](@cite) - - simple (integrated) sediment model described by [Soetaert2000](@citet), where organic matter - that sinks to the bottom is stored, decays into NO₃ and NH₄, and takes up O₂ in the process. - - Extended to attribute the corresponding release of DIC. - - Forced by O₂, NO₃, NH₄ and particle concentration in bottom cell. - - Instant remineralisation - - simple model from [Aumont2015](@citet), where sinking organic matter is instantly remineralised - and returned to the bottom cell - - some fraction is stored permanently in the sediment at an efficiency given by [RemineralisationFraction](@citet) -""" -module Boundaries - -export Sediments, GasExchange - -using Roots, Oceananigans, KernelAbstractions -using Oceananigans.Units -using Oceananigans.Architectures: device -using Oceananigans.Operators: volume, Az -using Oceananigans.BoundaryConditions: FluxBoundaryCondition -using Oceananigans.Fields: Field, fractional_indices -using Oceananigans.Grids: znodes - -import Adapt: adapt_structure, adapt -import Base: show, summary - -include("gasexchange.jl") -include("Sediments/Sediments.jl") - -end #module diff --git a/src/Boundaries/gasexchange.jl b/src/Boundaries/gasexchange.jl deleted file mode 100644 index b3b713b64..000000000 --- a/src/Boundaries/gasexchange.jl +++ /dev/null @@ -1,257 +0,0 @@ -##### -##### Carbonate chemistry to determine pCO₂ -##### As per OCMIP Phase 2 http://ocmip5.ipsl.jussieu.fr/OCMIP/phase2/simulations/Abiotic/Cchem/co2calc.f simplified as in PISCESv2 -##### - -using Oceananigans.Fields: fractional_indices, indices - -struct pCO₂{P0, P1, P2, PB, PW, FT} - solubility :: P0 - bicarbonate_dissociation :: P1 - carbonate_dissociation :: P2 - boric_acid_dissociation :: PB - water_dissociaiton :: PW # don't think this is the right name for it - lower_pH_bound :: FT - upper_pH_bound :: FT - boron_ratio :: FT - thermal_expansion :: FT - haline_contraction :: FT - - - pCO₂(solubility::P0, - bicarbonate_dissociation::P1, - carbonate_dissociation::P2, - boric_acid_dissociation::PB, - water_dissociaiton::PW, - lower_pH_bound::FT, - upper_pH_bound::FT, - boron_ratio::FT, - thermal_expansion::FT, - haline_contraction::FT) where {P0, P1, P2, PB, PW, FT} = - new{P0, P1, P2, PB, PW, FT}(solubility, - bicarbonate_dissociation, - carbonate_dissociation, - boric_acid_dissociation, - water_dissociaiton, - lower_pH_bound, upper_pH_bound, - boron_ratio, - thermal_expansion, haline_contraction) -end - -adapt_structure(to, pCO₂_model::pCO₂) = pCO₂(adapt(to, pCO₂_model.solubility), - adapt(to, pCO₂_model.bicarbonate_dissociation), - adapt(to, pCO₂_model.carbonate_dissociation), - adapt(to, pCO₂_model.boric_acid_dissociation), - adapt(to, pCO₂_model.water_dissociaiton), - adapt(to, pCO₂_model.lower_pH_bound), - adapt(to, pCO₂_model.upper_pH_bound), - adapt(to, pCO₂_model.boron_ratio), - adapt(to, pCO₂_model.thermal_expansion), - adapt(to, pCO₂_model.haline_contraction)) - -@inline function titrate_alkalinity(H, p) - return p.DIC * (p.k¹ * H + 2 * p.k¹ * p.k²) / (H ^ 2 + p.k¹ * H + p.k¹ * p.k²) - - (p.Alk - p.kʷ / H - p.boron / (1 + H / p.kᵇ)) -end - -@inline function (p::pCO₂)(DIC, Alk, T, S) - ρₒ = 1027 * (1 - p.thermal_expansion * T + p.haline_contraction * S) - - T += 273.15 - - Alk *= 1e-3 / ρₒ - DIC *= 1e-3 / ρₒ - - pk⁰ = p.solubility - pk¹ = p.bicarbonate_dissociation - pk² = p.carbonate_dissociation - pkᵇ = p.boric_acid_dissociation - pkʷ = p.water_dissociaiton - - k¹ = 10 ^ (pk¹.C + pk¹.invT / T + pk¹.logT * log(T) + pk¹.S * S + pk¹.S² * S ^ 2) - - k² = 10 ^ (pk².C + pk².invT / T + pk².S * S + pk².S² * S ^ 2) - - kᵇ = exp(pkᵇ.C + (pkᵇ.invT + pkᵇ.invTsqrtS * sqrt(S) + pkᵇ.invTS * S + pkᵇ.invTS¹⁵ * S ^ 1.5 + pkᵇ.invTS² * S ^ 2) / T - + pkᵇ.sqrtS * sqrt(S) + pkᵇ.S * S + (pkᵇ.logT + pkᵇ.logTsqrtS * sqrt(S) + pkᵇ.logTS * S) * log(T) + pkᵇ.TsqrtS * sqrt(S) * T) - - kʷ = exp(pkʷ.C + pkʷ.invT / T + pkʷ.logT * log(T) + - (pkʷ.sqrtS + pkʷ.sqrtSinvT / T + pkʷ.sqrtSlogT * log(T)) * sqrt(S) - + pkʷ.S * S) - - boron = p.boron_ratio * S - - params = (; k¹, k², kᵇ, kʷ, DIC, Alk, boron) - - if titrate_alkalinity(10 ^ - p.upper_pH_bound, params) * titrate_alkalinity(10 ^ - p.lower_pH_bound, params) < 0 - H = find_zero(titrate_alkalinity, (10 ^ - p.upper_pH_bound, 10 ^ - p.lower_pH_bound), Bisection(); atol = 1e-10, p = params) - else - H = NaN - end - - ff = exp(pk⁰.C + pk⁰.invT / T + - pk⁰.ClogT * log(T * pk⁰.logCT) + pk⁰.T² * T ^ 2 + - S * (pk⁰.S + pk⁰.ST * T + pk⁰.ST² * T ^ 2)) - - CO₂ = DIC * H ^ 2/ (H ^ 2 + k¹ * H + k¹ * k²) - pCO₂ = (CO₂ / ff) * 10 ^ 6 - - return pCO₂ # μatm -end - -OCMIP_solubility = (C = -162.8301, invT = 218.2968 * 100, logCT = 1 / 100, ClogT = 90.9241, T² = - 1.47696 / (100 ^ 2), ST² = 0.0049867 / (100 ^ 2), ST = -0.025225 / 100, S = .025695)# - -OCMIP_bicarbonate_dissociation = (C = 62.008, S = 0.0118, S² = -0.000116, invT = -3670.7, logT = -9.7944) - -OCMIP_carbonate_dissociation = (C = -4.777, S = 0.0184, S² = -0.000118, invT = -1394.7, logT = 0.0) - -OCMIP_boric_acid_dissociation = (C = 148.0248, invT = -8966.9, invTsqrtS = -2890.53, invTS = -77.942, invTS¹⁵ = 1.728, invTS² = - 0.0996, - sqrtS = 137.1942, S = 1.62142, logT = - 24.4344, logTsqrtS = - 25.085, logTS = - 0.2474, TsqrtS = 0.053105) - -OCMIP_water_dissociaiton = (C = 148.9652, invT = - 13847.26, logT = - 23.6521, sqrtSinvT = 118.67, sqrtS = -5.977, sqrtSlogT = 1.0495, S = -0.01615) - -OCMIP_default = pCO₂(OCMIP_solubility, - OCMIP_bicarbonate_dissociation, - OCMIP_carbonate_dissociation, - OCMIP_boric_acid_dissociation, - OCMIP_water_dissociaiton, - 0.0, 14.0, - 0.000232 / 1.80655 / 10.811, - 1.67e-4, 7.80e-4) - -##### -##### Gas exchange model of [Wanninkhof1992](@citet) -##### -# TODO: Implement Ho et al. 2006 wind speed dependence - -k(T, uₐᵥ, Sc_params) = 0.39 * (0.01 / 3600) * uₐᵥ ^ 2 * (Sc(T, Sc_params) / 660) ^ (-0.5)# m/s, may want to add variable wind speed instead of average wind here at some point -Sc(T, params) = params.A - params.B * T + params.C * T ^ 2 - params.D * T ^ 3 - -α(T, S, β_params) = β(T + 273.15, S, β_params) * (T + 273.15) * 0.00367#/(T+273.15) - disagree with origional paper but this matches dimensionless Henry coefficiet of 3.2x10⁻² at 298.15K, S=0. See https://www.wikiwand.com/en/Henry%27s_law -β(T, S, params) = exp(params.A₁ + params.A₂ * (100 / T) + params.A₃ * log(T / 100) + S * (params.B₁ + params.B₂ * (T / 100) + params.B₃ * (T / 100) ^ 2)) - -#now fairly sure that this is giving the correct result for CO₂ as βρ is the henrys coefficient which sould be ∼34mol/m³ atm -K(T, S, uₐᵥ, Sc_params, β_params, ρₒ) = k(T, uₐᵥ, Sc_params) * β(T + 273.15, S, β_params) * ρₒ #L=ρ\_wK₀ -> https://reader.elsevier.com/reader/sd/pii/0304420374900152 and here K₀=β - -##### -##### Boundary condition setup -##### - -struct GasExchange{G, ScP, βP, FT, AC, AP, PCO} <: Function - gas :: G - - schmidt_params :: ScP - solubility_params :: βP - ocean_density :: FT - air_concentration :: AC - air_pressure :: AP - average_wind_speed :: FT - - pCO₂ :: PCO -end - -adapt_structure(to, gasexchange::GasExchange) = GasExchange(adapt(to, gasexchange.gas), - adapt(to, gasexchange.schmidt_params), - adapt(to, gasexchange.solubility_params), - gasexchange.ocean_density, - adapt(to, gasexchange.air_concentration), - adapt(to, gasexchange.air_pressure), - gasexchange.average_wind_speed, - adapt(to, gasexchange.pCO₂)) - -""" - GasExchange(; gas, - schmidt_params::ScP = (CO₂ = (A=2073.1, B=125.62, C=3.6276, D=0.043219), - O₂ = (A=1953.4, B=128.0, C=3.9918, D=0.050091))[gas], - solubility_params::βP = (CO₂ = (A₁=-60.2409, A₂=93.4517, A₃=23.3585, B₁=0.023517, B₂=-0.023656, B₃=0.0047036), - O₂ = (A₁=-58.3877, A₂=85.8079, A₃=23.8439, B₁=-0.034892, B₂=0.015568, B₃=-0.0019387))[gas], - ocean_density::FT = 1026, # kg/m³ - air_concentration::AC = (CO₂ = 413.4, O₂ = 9352.7)[gas], # ppmv, mmolO₂/m³ (20.95 mol O₂/mol air, 0.0224m^3/mol air) - air_pressure::FT = 1.0, # atm - average_wind_speed::FT = 10, # m/s - field_dependencies = (CO₂ = (:DIC, :ALK), O₂ = (:OXY, ))[gas]) - -Construct an Oceananigans `FluxBoundaryCondition` for the exchange of `gas` with the relevant tracer (i.e., DIC for CO₂ and oxygen for O₂). -Please see note for other gases. - -Keyword arguments -================= - -- `gas`: (required) the gas to be exchanged, if `:CO₂` or `:O₂` are specified then all other settings may be infered -- `schmidt_params` : named tuple of parameters for calculating the Schmidt number using the parameterisation of [Wanninkhof1992](@citet) -- `solubility_params` : named tuple of parameters for calculating the solubility (for O₂ the Bunsen solubility and CO₂ K₀, see note) -- `ocean_density` : density of the ocean in kg/m³ -- `air_concentratio` : concentration of the gas in air in relivant units (i.e. ppmv for CO₂ and mmol O₂/m³ for O₂), can also be a function of x, y, t, or a field -- `air_pressure` : air pressure in atm (only used for CO₂), can also be a function of x, y, t, or a field -- `average_wind_speed` : average wind speed at 10m used to calculate the gas transfer velocity by the [Wanninkhof1992](@citet) parameterisation -- `field_dependencies` : tracer fields that gas exchange depends on, if the defaults have different names in your model you can specify as long as they are in the same order -- `pCO₂` : pCO₂ calculator - -!!! note "Gases other than CO₂ and O₂" - This model is fully capable of exchanging any gas but the parameters have only been configured for CO₂ and O₂, and the specific formulation - is only ensured for these gasses. For any gas where the [Wanninkhof1992](@citet) parameterisation returns the Bunsen Solubility Coefficient - this model will work out of the box and can just be passed new parameters. For the other solubility types (i.e. K₀, K' and f) you will need - to overload the `(gasexchange::GasExchange)` function to ensure the correct formulaiton. -""" -function GasExchange(; gas, - schmidt_params::ScP = (CO₂ = (A = 2073.1, B = 125.62, C = 3.6276, D = 0.043219), - O₂ = (A = 1953.4, B = 128.0, C = 3.9918, D = 0.050091))[gas], - solubility_params::βP = (CO₂ = (A₁ = -60.2409, A₂ = 93.4517, A₃ = 23.3585, B₁ = 0.023517, B₂ = -0.023656, B₃ = 0.0047036), - O₂ = (A₁ = -58.3877, A₂ = 85.8079, A₃ = 23.8439, B₁ = -0.034892, B₂ = 0.015568, B₃ = -0.0019387))[gas], - ocean_density::FT = 1024.5, # kg/m³ - air_concentration::AC = (CO₂ = 413.4, O₂ = 9352.7)[gas], # ppmv, mmolO₂/m³ (20.95 mol O₂/mol air, 0.0224m^3/mol air) - air_pressure::AP = 1.0, # atm - average_wind_speed::FT = 10.0, # m/s - field_dependencies = (CO₂ = (:DIC, :Alk), O₂ = (:O₂, ))[gas], - pCO₂::PCO = gas == :CO₂ ? OCMIP_default : nothing) where {ScP, βP, FT, AC, AP, PCO} - - gas = Val(gas) - G = typeof(gas) - - gasexchange = GasExchange{G, ScP, βP, FT, AC, AP, PCO}(gas, - schmidt_params, - solubility_params, - ocean_density, - air_concentration, - air_pressure, - average_wind_speed, - pCO₂) - - return FluxBoundaryCondition(gasexchange, field_dependencies = (field_dependencies..., :T, :S)) -end - -@inline function (gasexchange::GasExchange)(x, y, t, DIC, ALK, T, S) - conc = gasexchange.pCO₂(DIC, ALK, T, S) - return gasexchange(x, y, t, conc, T, S) -end - -@inline function (gasexchange::GasExchange)(x, y, t, conc, T, S) - return k(T, gasexchange.average_wind_speed, gasexchange.schmidt_params) * (conc - α(T, S, gasexchange.solubility_params) * get_value(x, y, t, gasexchange.air_concentration)) -end - -@inline function (gasexchange::GasExchange{<:Val{:CO₂}, <:Any, <:Any, <:Any, <:Any, <:Any})(x, y, t, conc, T, S) - return K(T, S, gasexchange.average_wind_speed, gasexchange.schmidt_params, gasexchange.solubility_params, gasexchange.ocean_density) * (conc - get_value(x, y, t, gasexchange.air_concentration) * get_value(x, y, t, gasexchange.air_pressure)) / 1000#μmol/m²s to mmolC/m²s not sure this is correct -end - -@inline get_value(x, y, t, air_concentration::Number) = air_concentration -@inline get_value(x, y, t, air_concentration::Function) = air_concentration(x, y, t) -#=It is not possible for this to work on GPU since we need the grid and locs before (since fields are reduced to vectors) - We probably need a more generic and better way todo this but here is not the place. - For now if you wanted to use data you could do similar to https://github.com/OceanBioME/GlobalOceanBioME.jl/blob/main/src/one_degree_surface_par.jl -# interpolate doesn't really work on 2D fields -@inline function get_value(x, y, t, conc::Field{LX, LY, LZ}) where {LX, LY, LZ} - grid = conc.grid - - i, j, _ = fractional_indices((x, y, 0.0), grid, LX(), LY(), Center()) - - ξ, i = mod(i, 1), Base.unsafe_trunc(Int, i) - η, j = mod(j, 1), Base.unsafe_trunc(Int, j) - - _, _, ks = indices(conc) - - return @inbounds ((1 - ξ) * (1 - η) * conc[i, j, ks[1]] - + (1 - ξ) * η * conc[i, j+1, ks[1]] - + ξ * (1 - η) * conc[i+1, j, ks[1]] - + ξ * η * conc[i+1, j+1, ks[1]]) -end -=# \ No newline at end of file diff --git a/src/BoxModel/boxmodel.jl b/src/BoxModel/boxmodel.jl index 48abc5c8f..f04800dd0 100644 --- a/src/BoxModel/boxmodel.jl +++ b/src/BoxModel/boxmodel.jl @@ -46,7 +46,7 @@ end forcing = NamedTuple(), timestepper = :RungeKutta3, clock = Clock(; time = 0.0), - prescribed_tracers = (:T, )) + prescribed_tracers::PT = (T = (t) -> 0, )) Constructs a box model of a `biogeochemistry` model. Once this has been constructed you can set initial condiitons by `set!(model, X=1.0...)`. @@ -57,7 +57,7 @@ Keyword Arguments - `forcing`: NamedTuple of additional forcing functions for the biogeochemical tracers to be integrated - `timestepper`: Timestepper to integrate model - `clock`: Oceananigans clock to keep track of time -- `prescribed_tracers`: Tuple of fields names (Symbols) which are not integrated but provided in `forcing` as a function of time with signature `f(t)` +- `prescribed_tracers`: named tuple of tracer names and function (`f(t)`) prescribing tracer values """ function BoxModel(; biogeochemistry::B, grid = BoxModelGrid(), @@ -84,9 +84,9 @@ end function update_state!(model::BoxModel, callbacks=[]; compute_tendencies = true) t = model.clock.time - for field in model.prescribed_tracers - if field in keys(model.fields) - @inbounds model.fields[field][1, 1, 1] = @inbounds model.forcing[field](t) + for (name, forcing) in pairs(model.prescribed_tracers) + if name in keys(model.fields) + @inbounds model.fields[name][1, 1, 1] = @inbounds forcing(t) end end diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index b24017e8f..89fa533db 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -49,7 +49,7 @@ using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers using OceanBioME.BoxModels: BoxModel -using OceanBioME.Boundaries.Sediments: sinking_flux +using OceanBioME.Models.Sediments: sinking_flux using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry @@ -64,7 +64,7 @@ import OceanBioME: maximum_sinking_velocity import Adapt: adapt_structure, adapt import Base: show, summary -import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers +import OceanBioME.Models.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry phytoplankton_preference :: FT diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 1104b1177..03139e24d 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -25,7 +25,7 @@ using Oceananigans.Fields: ZeroField using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR using OceanBioME: setup_velocity_fields, show_sinking_velocities using OceanBioME.BoxModels: BoxModel -using OceanBioME.Boundaries.Sediments: sinking_flux +using OceanBioME.Models.Sediments: sinking_flux import OceanBioME: redfield, conserved_tracers import Base: show, summary @@ -36,7 +36,7 @@ import Oceananigans.Biogeochemistry: required_biogeochemical_tracers, import OceanBioME: maximum_sinking_velocity -import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers +import OceanBioME.Models.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers import Adapt: adapt_structure, adapt diff --git a/src/Models/CarbonChemistry/CarbonChemistry.jl b/src/Models/CarbonChemistry/CarbonChemistry.jl new file mode 100644 index 000000000..63b6f11f8 --- /dev/null +++ b/src/Models/CarbonChemistry/CarbonChemistry.jl @@ -0,0 +1,14 @@ +""" +`CarbonChemistryModel` to solve chemical equilibrium parameterisations +""" +module CarbonChemistryModel + +export CarbonChemistry + +import Base: show, summary + +include("equilibrium_constants.jl") +include("carbon_chemistry.jl") +include("calcite_concentration.jl") + +end \ No newline at end of file diff --git a/src/Models/CarbonChemistry/calcite_concentration.jl b/src/Models/CarbonChemistry/calcite_concentration.jl new file mode 100644 index 000000000..9fca773f2 --- /dev/null +++ b/src/Models/CarbonChemistry/calcite_concentration.jl @@ -0,0 +1,81 @@ +function carbonate_concentration(cc::CarbonChemistry; + DIC, T, S, Alk = 0, pH = nothing, + P = nothing, + lon = 0, + lat = 0, + boron = 0.000232 / 10.811 * S / 1.80655, + sulfate = 0.14 / 96.06 * S / 1.80655, + fluoride = 0.000067 / 18.9984 * S / 1.80655, + silicate = 0, + phosphate = 0, + upper_pH_bound = 14, + lower_pH_bound = 0) + + ρₒ = cc.density_function(T, S, ifelse(isnothing(P), 0, P), lon, lat) + + # Centigrade to kelvin + T += 273.15 + + # mili-equivalents / m³ to equivalents / kg + Alk *= 1e-3 / ρₒ + + # mmol / m³ to mol / kg + DIC *= 1e-3 / ρₒ + phosphate *= 1e-3 / ρₒ + silicate *= 1e-3 / ρₒ + + # ionic strength + Is = cc.ionic_strength(S) + + # compute equilibrium constants + K1 = cc.carbonic_acid.K1(T, S; P) + K2 = cc.carbonic_acid.K2(T, S; P) + KB = cc.boric_acid(T, S; P) + KW = cc.water(T, S; P) + KS = cc.sulfate(T, S, Is; P) + KF = cc.fluoride(T, S, Is, KS; P) + KP1 = cc.phosphoric_acid.KP1(T, S; P) + KP2 = cc.phosphoric_acid.KP2(T, S; P) + KP3 = cc.phosphoric_acid.KP3(T, S; P) + KSi = cc.silicic_acid(T, S, Is) + + params = (; DIC, Alk, boron, sulfate, fluoride, silicate, phosphate, + K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi) + + # solve equilibrium for hydrogen ion concentration + H = solve_for_H(pH, params, upper_pH_bound, lower_pH_bound) + + # compute the calcite concentration + denom1 = (H * (H + K1)) + denom2 = (1 + K1 * K2 / denom1) + + return DIC * K1 * K2 / denom1 / denom2 +end + +function carbonate_saturation(cc::CarbonChemistry; + DIC, T, S, Alk = 0, pH = nothing, + P = nothing, + boron = 0.000232 / 10.811 * S / 1.80655, + sulfate = 0.14 / 96.06 * S / 1.80655, + fluoride = 0.000067 / 18.9984 * S / 1.80655, + calcium_ion_concentration = 0.0103 * S / 35, + silicate = 0, + phosphate = 0, + upper_pH_bound = 14, + lower_pH_bound = 0) + + CO₃²⁻ = carbonate_concentration(cc; + DIC, Alk, T, S, pH, + P, + boron, + sulfate, + fluoride, + silicate, + phosphate, + upper_pH_bound, + lower_pH_bound) + + KSP = cc.calcite_solubility(T, S; P) + + return calcium_ion_concentration * CO₃²⁻ / KSP +end \ No newline at end of file diff --git a/src/Models/CarbonChemistry/carbon_chemistry.jl b/src/Models/CarbonChemistry/carbon_chemistry.jl new file mode 100644 index 000000000..48eeeedf9 --- /dev/null +++ b/src/Models/CarbonChemistry/carbon_chemistry.jl @@ -0,0 +1,191 @@ +using Roots +using OceanBioME.Models: teos10_polynomial_approximation + +""" + CarbonChemistry(; ionic_strength = IonicStrength(), + solubility = K0(), + carbonic_acid = (K1 = K1(), K2 = K2()), + boric_acid = KB(), + water = KW(), + sulfate = KS(; ionic_strength), + fluoride = KF(; ionic_strength), + phosphoric_acid = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()), + silicic_acid = KSi(; ionic_strength), + calcite_solubility = KSP_calcite(), + density_function = teos10_polynomial_approximation) + +Carbon chemistry model capable of solving for sea water pCO₂ from DIC and +total alkalinity or DIC and pH. + +Default form from Dickson, A.G., Sabine, C.L. and Christian, J.R. (2007), +Guide to Best Practices for Ocean CO 2 Measurements. PICES Special Publication 3, 191 pp. + +See each parameters documentation for origional sources. + +Example +======= + +```jldoctest +julia> using OceanBioME + +julia> carbon_chemistry = CarbonChemistry() +`CarbonChemistry` model which solves for pCO₂ and pH + +julia> pCO₂ = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35) +1308.0843992121615 + +julia> pH = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35, return_pH = true) +7.502534641304366 + +julia> pCO₂_higher_pH = carbon_chemistry(; DIC = 2000, T = 10, S = 35, pH = 7.5) +1315.6558976217746 + +``` +""" +@kwdef struct CarbonChemistry{P0, PC, PB, PS, PF, PP, PSi, PW, IS, PKS, PRho} + ionic_strength :: IS = IonicStrength() + solubility :: P0 = K0() + carbonic_acid :: PC = (K1 = K1(), K2 = K2()) + boric_acid :: PB = KB() + water :: PW = KW() + sulfate :: PS = KS(; ionic_strength) + fluoride :: PF = KF(; ionic_strength, sulfate_constant = sulfate) + phosphoric_acid :: PP = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()) + silicic_acid :: PSi = KSi(; ionic_strength) + calcite_solubility :: PKS = KSP_calcite() + density_function :: PRho = teos10_polynomial_approximation +end + +""" + alkalinity_residual(H, p) + +Returns the difference between total alkalinity computed from `H`` (hydrogen ion +concentration), `DIC`, `borate`, `sulfate`, `phosphate`, `silicate`, and `fluoride` +concentration and chemical equilibrium constants specified in `p`, and the specified +total `Alk`alinity. + + TAlk = [HCO₃⁻] + 2[CO₃²⁻] + [B(OH)₄⁻] + [OH⁻] + [HPO₄²⁻] + 2[PO₄³⁻] + [SiO(OH)₃⁻] + + [NH₃] + [HS⁻] - [H⁺] - [HSO₄⁻] - [HF] - [H₃PO₄] + minor acids and bases + +Concentrations diagnosed as specified in Dickson et. al best practice descried in +`CarbonChemistry` docstring. + +Note ammonia (NH₃) is not currently included. +""" +@inline function alkalinity_residual(H, p) + carbonate_denom = H^2 + p.K1 * H + p.K1 * p.K2 + phosphorus_denom = H^3 + p.KP1 * H^2 + p.KP1 * p.KP2 * H + p.KP1 * p.KP2 * p.KP3 + sulfate_denom = 1 + p.sulfate / p.KS + + bicarbonate = p.K1 * H * p.DIC / carbonate_denom + carbonate = 2 * p.DIC * p.K1 * p.K2 / carbonate_denom + borate = p.boron / (1 + H / p.KB) + hydroxide = p.KW / H + hydrogen_phosphate = p.phosphate * p.KP1 * p.KP2 * H / phosphorus_denom + phosphate = 2 * p.phosphate * p.KP1 * p.KP2 * p.KP3 / phosphorus_denom + silicate = p.silicate / (1 + H / p.KSi) + free_hydrogen = - H / sulfate_denom + hydrogen_suplfate = - p.sulfate / (1 + p.KS / H / sulfate_denom) + hydrogen_fluoride = -p.fluoride / (1 + p.KF / H) + phosphoric_acid = -p.phosphate * H^3 / phosphorus_denom + + return (bicarbonate + + carbonate + + borate + + hydroxide + + hydrogen_phosphate + + phosphate + + silicate + + free_hydrogen + + hydrogen_suplfate + + hydrogen_fluoride + + phosphoric_acid + - p.Alk) +end + +""" + (p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing, + return_pH = false, + boron = 0.000232 / 10.811 * S / 1.80655, + sulfate = 0.14 / 96.06 * S / 1.80655, + fluoride = 0.000067 / 18.9984 * S / 1.80655, + silicate = 0, + phosphate = 0, + upper_pH_bound = 14, + lower_pH_bound = 0) + +Calculates `fCO₂` in sea water with `DIC`, `Alk`alinity, `T`emperature, and `S`alinity +unless `pH` is specified, in which case intermediate computation of `pH` is skipped and +`pCO₂` is calculated from the `DIC`, `T`, `S` and `pH`. + +Alternativly, `pH` is returned if `return_pH` is `true`. +""" +@inline function (p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing, + P = nothing, + lon = 0, + lat = 0, + return_pH = false, + boron = 0.000232 / 10.811 * S / 1.80655, + sulfate = 0.14 / 96.06 * S / 1.80655, + fluoride = 0.000067 / 18.9984 * S / 1.80655, + silicate = 0, + phosphate = 0, + upper_pH_bound = 14, + lower_pH_bound = 0) + + ρₒ = p.density_function(T, S, ifelse(isnothing(P), 0, P), lon, lat) + + # Centigrade to kelvin + T += 273.15 + + # mili-equivalents / m³ to equivalents / kg + Alk *= 1e-3 / ρₒ + + # mmol / m³ to mol / kg + DIC *= 1e-3 / ρₒ + phosphate *= 1e-3 / ρₒ + silicate *= 1e-3 / ρₒ + + # ionic strength + Is = p.ionic_strength(S) + + # compute equilibrium constants + K1 = p.carbonic_acid.K1(T, S; P) + K2 = p.carbonic_acid.K2(T, S; P) + KB = p.boric_acid(T, S; P) + KW = p.water(T, S; P) + KS = p.sulfate(T, S, Is; P) + KF = p.fluoride(T, S, Is, KS; P) + KP1 = p.phosphoric_acid.KP1(T, S; P) + KP2 = p.phosphoric_acid.KP2(T, S; P) + KP3 = p.phosphoric_acid.KP3(T, S; P) + KSi = p.silicic_acid(T, S, Is) + + params = (; DIC, Alk, boron, sulfate, fluoride, silicate, phosphate, + K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi) + + # solve equilibrium for hydrogen ion concentration + H = solve_for_H(pH, params, upper_pH_bound, lower_pH_bound) + + # compute solubility equilibrium constant + K0 = p.solubility(T, S) + + # compute pCO₂ + CO₂ = DIC * H ^ 2 / (H ^ 2 + K1 * H + K1 * K2) + fCO₂ = (CO₂ / K0) * 10 ^ 6 # μatm + + # compute pH + pH = -log10(H) + + return ifelse(return_pH, pH, fCO₂) +end + +# solves `alkalinity_residual` for pH +solve_for_H(pH, args...) = 10.0 ^ - pH + +solve_for_H(::Nothing, params, upper_pH_bound, lower_pH_bound) = + find_zero(alkalinity_residual, (10.0 ^ - upper_pH_bound, 10.0 ^ - lower_pH_bound), Bisection(); atol = 1e-10, p = params) + +# display +summary(::IO, ::CarbonChemistry) = string("`CarbonChemistry` model") +show(io::IO, ::CarbonChemistry) = print(io, "`CarbonChemistry` model which solves for pCO₂ and pH") \ No newline at end of file diff --git a/src/Models/CarbonChemistry/equilibrium_constants.jl b/src/Models/CarbonChemistry/equilibrium_constants.jl new file mode 100644 index 000000000..07e3856b7 --- /dev/null +++ b/src/Models/CarbonChemistry/equilibrium_constants.jl @@ -0,0 +1,721 @@ +""" + PressureCorrection(a₀, a₁, a₂, + b₀, b₁, b₂, + R) + +Parameterisation for the pressure effect on thermodynamic constants. + +Form from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). +""" +@kwdef struct PressureCorrection{FT} + a₀ :: FT + a₁ :: FT + a₂ :: FT + b₀ :: FT + b₁ :: FT + + R :: FT = 83.14472 +end + +@inline function (pc::PressureCorrection)(Tk, P) + Tc = Tk - 273.15 + + ΔV = pc.a₀ + pc.a₁ * Tc + pc.a₂ * Tc^2 + Δκ = pc.b₀ + pc.b₁ * Tc + + RT = pc.R * Tk + + return exp((-ΔV + 0.5 * Δκ * P) * P / RT) +end + +@inline (pc::PressureCorrection)(T, ::Nothing) = 1 + +summary(::IO, ::PressureCorrection) = string("Equilibrium constant pressure correction") +show(io::IO, pc::PressureCorrection) = print(io, "Equilibrium constant pressure correction\n", + " ΔV = $(pc.a₀) + $(pc.a₁) Tc + $(pc.a₂) * Tc²\n", + " Δκ = $(pc.b₀) + $(pc.b₁) T\n", + " ln(kᵖ/k⁰) = (-ΔV + Δκ P / 2) * P / RT") + +""" + K0(; constant = -60.2409, + inverse_T = 93.4517 * 100, + log_T = 23.3585, + T² = 0.0, + S = 0.023517, + ST = -0.023656 / 100, + ST² = 0.0047036 / 100^2) + +Parameterisation for carbon dioxide solubility equilibrium constant. + + CO₂(g) ⇌ CO₂*(aq) + + K₀ = [CO₂*(aq)]/f(CO₂) + +Default values from Weiss, R.F. (1974, Mar. Chem., 2, 203–215). +""" +@kwdef struct K0{FT} + constant :: FT = -60.2409 + inverse_T :: FT = 93.4517 * 100 + log_T :: FT = 23.3585 + T² :: FT = 0.0 + S :: FT = 0.023517 + ST :: FT = -0.023656 / 100 + ST² :: FT = 0.0047036 / 100^2 +end + +@inline (c::K0)(T, S; P = nothing) = + exp(c.constant + + c.inverse_T / T + + c.log_T * (log(T) - log(100)) + + c.T² * T^2 + + (c.S + c.ST * T + c.ST² * T^2) * S) + +summary(::IO, ::K0) = string("Solubility constant") +show(io::IO, k0::K0) = print(io, "Solubility constant\n", + " ln(k₀/k°) = $(k0.constant) + $(k0.inverse_T) / T + $(k0.log_T) (log(T) - log(100)) + $(k0.T²) T² + ($(k0.S) + $(k0.ST) T + $(k0.ST²) T²)S") + +""" + K1(; constant = 61.2172, + inverse_T = -3633.86, + log_T = -9.67770, + S = 0.011555, + S² = -0.0001152, + pressure_correction = + PressureCorrection(; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877)) + +Parameterisation for aquious carbon dioxide - bicarbonate dissociation equilibrium constant. + + CO₂*(aq) + H₂O ⇌ H₂CO₃ ⇌ HCO₃⁻ + H⁺ + + K₁ = [H⁺][HCO₃⁻]/[CO₂*] + +Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119). +""" +@kwdef struct K1{FT, PC} + constant :: FT = 61.2172 + inverse_T :: FT = -3633.86 + log_T :: FT = -9.67770 + S :: FT = 0.011555 + S² :: FT = -0.0001152 + + pressure_correction :: PC = + PressureCorrection(; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877) +end + +@inline (c::K1)(T, S; P = nothing) = + c.pressure_correction(T, P) * + 10 ^ (c.constant + c.inverse_T / T + c.log_T * log(T) + c.S * S + c.S² * S^2) + +summary(::IO, ::K1) = string("First carbon dioxide dissociation constant") +show(io::IO, k1::K1) = print(io, "First carbon dioxide dissociation constant\n", + " log₁₀(k₁/k°) = $(k1.constant) + $(k1.inverse_T) / T + $(k1.log_T) log(T) + $(k1.S) S + $(k1.S²) S²") + +""" + K2(; constant = -25.9290, + inverse_T = -471.78, + log_T = 3.16967, + S = 0.01781, + S² = -0.0001122, + pressure_correction = + PressureCorrection(; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475)) + +Parameterisation for bicarbonate dissociation equilibrium constant. + + HCO₃⁻ ⇌ CO₃²⁻ + H⁺ + + K₂ = [H⁺][CO₃²⁻]/[HCO₃⁻] + +Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119). +""" +@kwdef struct K2{FT, PC} + constant :: FT = -25.9290 + inverse_T :: FT = -471.78 + log_T :: FT = 3.16967 + S :: FT = 0.01781 + S² :: FT = -0.0001122 + + pressure_correction :: PC = + PressureCorrection(; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475) +end + +@inline (c::K2)(T, S; P = nothing) = + c.pressure_correction(T, P) * + 10 ^ (c.constant + c.inverse_T / T + c.S * S + c.S² * S^2 + c.log_T * log(T)) + +summary(::IO, ::K2) = string("Second carbon dioxide dissociation constant") +show(io::IO, k2::K2) = print(io, "Second carbon dioxide dissociation constant\n", + " log₁₀(k₂/k°) = $(k2.constant) + $(k2.inverse_T) / T + $(k2.log_T) log(T) + $(k2.S) S + $(k2.S²) S²") + +""" + KB(; constant = 148.0248, + inverse_T = -8966.90, + inverse_T_sqrt_S = -2890.53, + inverse_T_S = -77.942, + inverse_T_sqrt_S³ = 1.728, + inverse_T_S² = -0.0996, + sqrt_S = 137.1942, + S = 1.62142, + log_T = -24.4344, + log_T_sqrt_S = -25.085, + S_log_T = -0.2474, + T_sqrt_S = 0.053105, + pressure_correction = + PressureCorrection(; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0)) + +Parameterisation for boric acid equilibrium with water. + + B(OH)₃ + H₂O ⇌ B(OH)₄⁻ + H⁺ + + Kᵇ = [H⁺][B(OH)₄⁻]/[B(OH)₃] + +Default values from Dickson (1990, Deep-Sea Res., 37, 755–766). +""" +@kwdef struct KB{FT, PC} + constant :: FT = 148.0248 + inverse_T :: FT = -8966.90 + inverse_T_sqrt_S :: FT = -2890.53 + inverse_T_S :: FT = -77.942 + inverse_T_sqrt_S³ :: FT = 1.728 + inverse_T_S² :: FT = -0.0996 + sqrt_S :: FT = 137.1942 + S :: FT = 1.62142 + log_T :: FT = -24.4344 + log_T_sqrt_S :: FT = -25.085 + S_log_T :: FT = -0.2474 + T_sqrt_S :: FT = 0.053105 + + pressure_correction :: PC = + PressureCorrection(; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0) +end + +@inline (c::KB)(T, S; P = nothing) = + c.pressure_correction(T, P) * + exp(c.constant + + (c.inverse_T + c.inverse_T_sqrt_S * √S + c.inverse_T_S * S + c.inverse_T_sqrt_S³ * S^1.5 + c.inverse_T_S² * S^2) / T + + c.sqrt_S * √S + + c.S * S + + (c.log_T + c.log_T_sqrt_S * √S + c.S_log_T * S ) * log(T) + + c.T_sqrt_S * √S * T) + +summary(::IO, ::KB) = string("Boric acid dissociation constant") +show(io::IO, c::KB) = print(io, "Boric acid dissociation constant\n", + " ln(kᵇ/k°) = $(c.constant) + ($(c.inverse_T) + $(c.inverse_T_sqrt_S) √S + $(c.inverse_T_S) S + $(c.inverse_T_sqrt_S³) √S³ + $(c.inverse_T_S²) S²) / T + + $(c.sqrt_S) * √S + + $(c.S) * S + + ($(c.log_T) + $(c.log_T_sqrt_S) √S + $(c.S_log_T) S ) * log(T) + + $(c.T_sqrt_S) * √S * T") + +""" + KW(; constant = 148.9652, + inverse_T = -13847.26, + log_T = -23.6521, + sqrt_S = -5.977, + inverse_T_sqrt_S = 118.67, + log_T_sqrt_S = 1.0495, + S = -0.01615, + pressure_correction = + PressureCorrection(; a₀=-20.02, a₁=0.1119, a₂=-0.001409, b₀=-0.00513, b₁=0.0000794)) + +Parameterisation for water dissociation equilibrium constant. + + H₂O ⇌ OH⁻ + H⁺ + + Kʷ = [H⁺][OH⁻] + +Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). +""" +@kwdef struct KW{FT, PC} + constant :: FT = 148.9652 + inverse_T :: FT = -13847.26 + log_T :: FT = -23.6521 + sqrt_S :: FT = -5.977 + inverse_T_sqrt_S :: FT = 118.67 + log_T_sqrt_S :: FT = 1.0495 + S :: FT = -0.01615 + + pressure_correction :: PC = + PressureCorrection(; a₀=-20.02, a₁=0.1119, a₂=-0.001409, b₀=-0.00513, b₁=0.0000794) +end + +@inline (c::KW)(T, S; P = nothing) = + c.pressure_correction(T, P) * + exp(c.constant + + c.inverse_T / T + + c.log_T * log(T) + + (c.sqrt_S + c.inverse_T_sqrt_S / T + c.log_T_sqrt_S * log(T))* √S + + c.S * S) + +summary(::IO, ::KW) = string("Water dissociation constant") +show(io::IO, c::KW) = print(io, "Water dissociation constant\n", + " ln(kʷ/k°) = $(c.constant) + + $(c.inverse_T) / T + + $(c.log_T) log(T) + + ($(c.sqrt_S) + $(c.inverse_T_sqrt_S) / T + $(c.log_T_sqrt_S) log(T)) √S + + $(c.S) * S") + + +""" + IonicStrength(; a = 19.924, + b = 1000.0, + c = -1.005) + +Parameterisation of the ionic strength of sea water. + + Is(S) = aS / (b + cS) + +Default values from Dickson (1990, Chem. Thermodyn., 22, 113–127). +""" +@kwdef struct IonicStrength{FT} + a :: FT = 19.924 + b :: FT = 1000.0 + c :: FT = -1.005 +end + +@inline (c::IonicStrength)(S) = c.a * S / (c.b + c.c * S) + +summary(::IO, ::IonicStrength) = string("Ionic strength") +show(io::IO, c::IonicStrength) = print(io, "Ionic strength\n", + " Is = $(c.a) S/($(c.b) + $(c.c) S)") + +""" + KS(; constant = 148.9652, + inverse_T = -13847.26, + log_T = -23.6521, + sqrt_S = -5.977, + inverse_T_sqrt_S = 118.67, + log_T_sqrt_S = 1.0495, + S = -0.01615, + pressure_correction = + PressureCorrection(; a₀=-18.03, a₁=0.0466, a₂=0.000316, b₀=-0.00453, b₁=0.00009)) + +Parameterisation for bisulfate dissociation equilibrium constant. + + HSO₄⁻ ⇌ SO₄²⁻ + H⁺ + + Kˢ = [H⁺][SO₄²⁻]/[HSO₄⁻] + +Default values from Dickson (1990, Chem. Thermodyn., 22, 113–127). +""" +@kwdef struct KS{IS, FT, PC} + ionic_strength :: IS = IonicStrength() + + constant :: FT = 141.328 + inverse_T :: FT = -4276.1 + log_T :: FT = -23.093 + sqrt_Is :: FT = 324.57 + inverse_T_sqrt_Is :: FT = -13856.0 + log_T_sqrt_Is :: FT = -47.986 + Is :: FT = -771.54 + inverse_T_Is :: FT = 35474.0 + log_T_Is :: FT = 114.723 + inverse_T_sqrt_Is³ :: FT = -2698.0 + inverse_T_Is² :: FT = 1776.0 + log_S :: FT = -0.001005 + + pressure_correction :: PC = + PressureCorrection(; a₀=-18.03, a₁=0.0466, a₂=0.000316, b₀=-0.00453, b₁=0.00009) +end + +@inline (c::KS)(T, S, Is = c.ionic_strength(S); P = nothing) = + c.pressure_correction(T, P) * + exp(c.constant + + c.inverse_T / T + + c.log_T * log(T) + + (c.sqrt_Is + c.inverse_T_sqrt_Is / T + c.log_T_sqrt_Is * log(T)) * √Is + + (c.Is + c.inverse_T_Is / T + c.log_T_Is * log(T)) * Is + + c.inverse_T_sqrt_Is³ * Is^1.5 / T + + c.inverse_T_Is² * Is^2 / T + + log(1 + c.log_S * S)) + +summary(::IO, ::KS) = string("Bisulfate dissociation constant") +show(io::IO, c::KS) = print(io, "Bisulfate dissociation constant\n", + " ln(kˢ/k°) = $(c.constant) + + $(c.inverse_T) / T + + $(c.log_T) log(T) + + ($(c.sqrt_Is) + $(c.inverse_T_sqrt_Is) / T + $(c.log_T_sqrt_Is) log(T)) √Is + + ($(c.Is) + $(c.inverse_T_Is) / T + $(c.log_T_Is) log(T)) Is + + $(c.inverse_T_sqrt_Is³) √Is³ / T + + $(c.inverse_T_Is²) Is² / T + + log(1 + $(c.log_S) S)") + +""" + KF(; ionic_strength = IonicStrength(), + sulfate_constant = KS(; ionic_strength), + constant = -9.68, + inverse_T = 874.0, + sqrt_S = 0.111, + log_S = 0.0, + log_S_KS = 0.0, + pressure_correction = + PressureCorrection(; a₀=-9.78, a₁=-0.0090, a₂=-0.000942, b₀=-0.00391, b₁=0.000054)) + +Parameterisation for hydrogen fluoride dissociation equilibrium constant. + + HF ⇌ F⁻ + H⁺ + + Kᶠ = [H⁺][F⁻]/[HF] + +Default values from Perez and Fraga (1987, Mar. Chem., 21, 161–168). +""" +@kwdef struct KF{IS, KS, FT, PC} + ionic_strength :: IS = IonicStrength() + sulfate_constant :: KS = KS(; ionic_strength) + + constant :: FT = -9.68 + inverse_T :: FT = 874.0 + sqrt_S :: FT = 0.111 + log_S :: FT = 0.0 + log_S_KS :: FT = 0.0 + + pressure_correction :: PC = + PressureCorrection(; a₀=-9.78, a₁=-0.0090, a₂=-0.000942, b₀=-0.00391, b₁=0.000054) +end + +@inline (c::KF)(T, S, Is = c.ionic_strength(S), KS = c.sulfate_constant(T, S, Is); P = nothing) = + c.pressure_correction(T, P) * + exp(c.constant + + c.inverse_T / T + + c.sqrt_S * √S + + log(1 + c.log_S * S) + + log(1 + c.log_S_KS * S / KS)) + +summary(::IO, ::KF) = string("Hydrogen fluoride dissociation constant") +show(io::IO, c::KF) = print(io, "Hydrogen fluoride dissociation constant\n", + " ln(kᶠ/k°) = $(c.constant) + + $(c.inverse_T) / T + + $(c.sqrt_S) √S + + log(1 + $(c.log_S) S) + + log(1 + $(c.log_S_KS) S / Kˢ)") + +""" + KP(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) + +Generic equilibrium constant parameterisation of the form used by +Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677) for phosphoric +acid dissociation. +""" +struct KP{FT, PC} + constant :: FT + inverse_T :: FT + log_T :: FT + sqrt_S :: FT + inverse_T_sqrt_S :: FT + S :: FT + inverse_T_S :: FT + + pressure_correction :: PC +end + +@inline (c::KP)(T, S; P = nothing) = + c.pressure_correction(T, P) * + exp(c.constant + + c.inverse_T / T + + c.log_T * log(T) + + (c.sqrt_S + c.inverse_T_sqrt_S / T) * √S + + (c.S + c.inverse_T_S / T) * S) + +summary(::IO, ::KP) = string("Phosphate dissociation constant") +show(io::IO, c::KP) = print(io, "Phosphate dissociation constant\n", + " ln(kᵖⁿ/k°) = $(c.constant) + + $(c.inverse_T) / T + + $(c.log_T) log(T) + + ($(c.sqrt_S) + $(c.inverse_T_sqrt_S) / T) √S + + ($(c.S) + $(c.inverse_T_S) / T) S") + +""" + KP1(; constant = 115.525, + inverse_T = -4576.752, + log_T = - 18.453, + sqrt_S = 0.69171, + inverse_T_sqrt_S = -106.736, + S = -0.01844, + inverse_T_S = -0.65643, + pressure_correction = + PressureCorrection(; a₀=-14.51, a₁=0.1211, a₂=-0.000321, b₀=-0.00267, b₁=0.0000427)) + +Instance of `KP` returning the first phosphocic acid equilibrium constant. + + H₃PO₄ ⇌ H₂PO₄⁻ + H⁺ + + Kᵖ¹ = [H⁺][H₂PO₄]/[H₃PO₄] + +Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). +""" +KP1(; constant = 115.525, + inverse_T = -4576.752, + log_T = - 18.453, + sqrt_S = 0.69171, + inverse_T_sqrt_S = -106.736, + S = -0.01844, + inverse_T_S = -0.65643, + pressure_correction = + PressureCorrection(; a₀=-14.51, a₁=0.1211, a₂=-0.000321, b₀=-0.00267, b₁=0.0000427)) = + KP(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) + +""" + KP2(; constant = 172.0883, + inverse_T = -8814.715, + log_T = -27.927, + sqrt_S = 1.3566, + inverse_T_sqrt_S = -160.340, + S = -0.05778, + inverse_T_S = 0.37335, + pressure_correction = + PressureCorrection(; a₀=-23.12, a₁=0.1758, a₂=-0.002647, b₀=-0.00515, b₁=0.00009)) + +Instance of `KP` returning the second phosphocic acid equilibrium constant. + + H₂PO₄⁻ ⇌ HPO₄²⁻ + H⁺ + + Kᵖ² = [H⁺][HPO₄²⁻]/[H₂PO₄⁻] + +Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). +""" +KP2(; constant = 172.0883, + inverse_T = -8814.715, + log_T = -27.927, + sqrt_S = 1.3566, + inverse_T_sqrt_S = -160.340, + S = -0.05778, + inverse_T_S = 0.37335, + pressure_correction = + PressureCorrection(; a₀=-23.12, a₁=0.1758, a₂=-0.002647, b₀=-0.00515, b₁=0.00009)) = + KP(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) + +""" + KP3(; constant = -18.141, + inverse_T = -3070.75, + log_T = 0.0, + sqrt_S = 2.81197, + inverse_T_sqrt_S = 17.27039, + S = -0.09984, + inverse_T_S = -44.99486, + pressure_correction = + PressureCorrection(; a₀=-26.57, a₁=0.2020, a₂=-0.0030420, b₀=-0.00408, b₁=0.0000714)) + +Instance of `KP` returning the third phosphocic acid equilibrium constant. + + HPO₄²⁻ ⇌ PO₄ + H⁺ + + Kᵖ³ = [H⁺][PO₄³⁻]/[HPO₄⁻] + +Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). +""" +KP3(; constant = - 18.141, + inverse_T = -3070.75, + log_T = 0.0, + sqrt_S = 2.81197, + inverse_T_sqrt_S = 17.27039, + S = -0.09984, + inverse_T_S = -44.99486, + pressure_correction = + PressureCorrection(; a₀=-26.57, a₁=0.2020, a₂=-0.0030420, b₀=-0.00408, b₁=0.0000714)) = + KP(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) + +""" + KSi(; ionic_strength = IonicStrength(), + constant = 117.385, + inverse_T = -8904.2, + log_T = -19.334, + sqrt_Is = 3.5913, + inverse_T_sqrt_Is = -458.79, + Is = -1.5998, + inverse_T_Is = 188.74, + Is² = 0.07871, + inverse_T_Is² = -12.1652, + log_S = -0.001005) + +Parameterisation for silicic acid dissociation equilibrium constant. + + Si(OH)₄ ⇌ SiO(OH)₃⁻ + H⁺ + + Kʷ = [H⁺][SiO(OH)₃⁻]/[Si(OH)₄] + +Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). +""" +@kwdef struct KSi{IS, FT} + ionic_strength :: IS = IonicStrength() + + constant :: FT = 117.385 + inverse_T :: FT = -8904.2 + log_T :: FT = -19.334 + sqrt_Is :: FT = 3.5913 + inverse_T_sqrt_Is :: FT = -458.79 + Is :: FT = -1.5998 + inverse_T_Is :: FT = 188.74 + Is² :: FT = 0.07871 + inverse_T_Is² :: FT = -12.1652 + log_S :: FT = -0.001005 +end + +@inline (c::KSi)(T, S, Is = c.ionic_strength(S); P = nothing) = + exp(c.constant + + c.inverse_T / T + + c.log_T * log(T) + + (c.sqrt_Is + c.inverse_T_sqrt_Is / T) * √Is + + (c.Is + c.inverse_T_Is / T) * Is + + (c.Is² + c.inverse_T_Is² / T) * Is^2 + + log(1 + c.log_S * S)) + +summary(::IO, ::KSi) = string("Silicic acid constant") +show(io::IO, c::KSi) = print(io, "Silicic acid constant\n", + " ln(kˢⁱ/k°) = $(c.constant) + + $(c.inverse_T) / T + + $(c.log_T) log(T) + + ($(c.sqrt_Is) + $(c.inverse_T_sqrt_Is) / T) √Is + + ($(c.Is) + $(c.inverse_T_Is) / T) Is + + ($(c.Is²) + $(c.inverse_T_Is²) / T) Is² + + log(1 + $(c.log_S) S)") + +""" + KSP(therm_constant, + therm_T, + therm_inverse_T, + therm_log_T, + sea_sqrt_S, + sea_T_sqrt_S, + sea_inverse_T_sqrt_S, + sea_S, + sea_S_sqrt_S³, + pressure_correction) + +Generic CaCO₃ solubility parameterisation of the form given by +Form from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). +""" +struct KSP{FT, PC} + therm_constant :: FT + therm_T :: FT + therm_inverse_T :: FT + therm_log_T :: FT + sea_sqrt_S :: FT + sea_T_sqrt_S :: FT + sea_inverse_T_sqrt_S :: FT + sea_S :: FT + sea_S_sqrt_S³ :: FT + + pressure_correction :: PC +end + +@inline function (c::KSP)(T, S; P = nothing) + + pressure_correction = c.pressure_correction(T, P) + + lnK_therm = c.therm_constant + c.therm_T * T + c.therm_inverse_T / T + c.therm_log_T * log10(T) # seems wrong + lnK_sea = ((c.sea_sqrt_S + c.sea_T_sqrt_S * T + c.sea_inverse_T_sqrt_S / T) * √S + + c.sea_S * S + + c.sea_S_sqrt_S³ * S^1.5) + + return pressure_correction * exp(lnK_therm + lnK_sea) +end + +summary(::IO, ::KSP) = string("Calcite solubility") +show(io::IO, c::KSP) = print(io, "Calcite solubility\n", + " ln(kₛₚ) = $(c.therm_constant) + $(c.therm_T) T + $(c.therm_inverse_T) / T + $(c.therm_log_T) log(T)\n", + " ln(kₛₚˢ) = ln(kₛₚ) + ($(c.sea_sqrt_S) + $(c.sea_T_sqrt_S) T + $(c.sea_inverse_T_sqrt_S) / T) √S\n", + " + $(c.sea_S) S + $(c.sea_S_sqrt_S³) √S³") + +""" + KSP_calcite(; therm_constant = -171.9065, + therm_T = -0.077993, + therm_inverse_T = 2839.319, + therm_log_T = 71.595, + sea_sqrt_S = -0.77712, + sea_T_sqrt_S = 0.0028426, + sea_inverse_T_sqrt_S = 178.34, + sea_S = -0.07711, + sea_S_sqrt_S³ = 0.0041249, + pressure_correction = + PressureCorrection(; a₀=-48.76, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) + +Instance of `KSP` returning calcite solubility. + +Default values from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). +""" +KSP_calcite(; therm_constant = -171.9065, + therm_T = -0.077993, + therm_inverse_T = 2839.319, + therm_log_T = 71.595, + sea_sqrt_S = -0.77712, + sea_T_sqrt_S = 0.0028426, + sea_inverse_T_sqrt_S = 178.34, + sea_S = -0.07711, + sea_S_sqrt_S³ = 0.0041249, + pressure_correction = + PressureCorrection(; a₀=-48.76, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = + KSP(therm_constant, + therm_T, + therm_inverse_T, + therm_log_T, + sea_sqrt_S, + sea_T_sqrt_S, + sea_inverse_T_sqrt_S, + sea_S, + sea_S_sqrt_S³, + pressure_correction) + +""" +KSP_aragonite(; therm_constant = -171.945, + therm_T = -0.077993, + therm_inverse_T = 2903.293, + therm_log_T = 71.595, + sea_sqrt_S = -0.068393, + sea_T_sqrt_S = 0.0017276, + sea_inverse_T_sqrt_S = 88.135, + sea_S = -0.10018, + sea_S_sqrt_S³ = 0.0059415, + pressure_correction = + PressureCorrection(; a₀=-45.96, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) + +Instance of `KSP` returning calcite solubility. + +Default values from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). +""" +KSP_aragonite(; therm_constant = -171.945, + therm_T = -0.077993, + therm_inverse_T = 2903.293, + therm_log_T = 71.595, + sea_sqrt_S = -0.068393, + sea_T_sqrt_S = 0.0017276, + sea_inverse_T_sqrt_S = 88.135, + sea_S = -0.10018, + sea_S_sqrt_S³ = 0.0059415, + pressure_correction = + PressureCorrection(; a₀=-45.96, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = + KSP(therm_constant, + therm_T, + therm_inverse_T, + therm_log_T, + sea_sqrt_S, + sea_T_sqrt_S, + sea_inverse_T_sqrt_S, + sea_S, + sea_S_sqrt_S³, + pressure_correction) \ No newline at end of file diff --git a/src/Models/GasExchange/GasExchange.jl b/src/Models/GasExchange/GasExchange.jl new file mode 100644 index 000000000..6df4ce011 --- /dev/null +++ b/src/Models/GasExchange/GasExchange.jl @@ -0,0 +1,138 @@ +""" +`GasExchangeModel` to solve chemical equilibrium parameterisations +""" +module GasExchangeModel + +export GasExchange, + CarbonDioxideGasExchangeBoundaryCondition, + OxygenGasExchangeBoundaryCondition, + GasExchangeBoundaryCondition, + ScaledGasTransferVelocity, + SchmidtScaledTransferVelocity, + CarbonDioxidePolynomialSchmidtNumber, + OxygenPolynomialSchmidtNumber + +using Adapt +using Oceananigans.BoundaryConditions: FluxBoundaryCondition +using Oceananigans.Fields: Center +using Oceananigans.Grids: xnode, ynode + +using OceanBioME.Models.CarbonChemistryModel: CarbonChemistry + +import Base: show, summary +import Adapt: adapt_structure + +const ATM = 101325 # Pa +const GAS_CONSTANT = 8.31446261815324 # J / kg / mol + +include("surface_values.jl") + +# just get the value from the model tracer +struct OxygenConcentration end + +summary(::OxygenConcentration) = "Model tracer `OxygenConcentration`" + +@inline surface_value(::OxygenConcentration, i, j, grid, clock, model_fields) = @inbounds model_fields.O₂[i, j, grid.Nz] + +# include all the bits +include("generic_parameterisations.jl") +include("gas_exchange.jl") +include("carbon_dioxide_concentration.jl") +include("schmidt_number.jl") +include("gas_transfer_velocity.jl") + +using .ScaledGasTransferVelocity + +# wrappers to produce boundary conditions + +""" + GasExchangeBoundaryCondition(; water_concentration, + air_concentration, + transfer_velocity, + wind_speed) + +Returns a `FluxBoundaryCondition` for the gas exchange between `water_concentration` and `air_concentration` +with `transfer_velocity`. + +`water_concentration`, `air_concentration` and `wind_speed` can either be numbers, +functions of the form `(x, y, t)`, functions of the form `(i, j, grid, clock, model_fields)` +if `discrete_form` is set to true, or any kind of `Field`. + +`water_concentration` should usually be a `[Tracer]Concentration` where is the name of the +tracer (you will have to build your own if this is not `OxygenConcentration`), +or a `CarbonDioxideConcentration` which diagnoses the partial pressure of CO₂ in the water. + +`transfer_velocity` should be a function of the form `k(u₁₀, T)`. +""" +function GasExchangeBoundaryCondition(; water_concentration, + air_concentration, + transfer_velocity, + wind_speed, + discrete_form = false) + + wind_speed = normalise_surface_function(wind_speed; discrete_form) + air_concentration = normalise_surface_function(air_concentration; discrete_form) + + exchange_function = GasExchange(wind_speed, transfer_velocity, water_concentration, air_concentration) + + return FluxBoundaryCondition(exchange_function; discrete_form = true) +end + +""" + CarbonDioxideGasExchangeBoundaryCondition(; carbon_chemistry = CarbonChemistry(), + transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = CarbonDioxidePolynomialSchmidtNumber()), + air_concentration = 413, # ppmv + wind_speed = 2, + water_concentration = nothing, + silicate_and_phosphate_names = nothing, + kwargs...) + +Returns a `FluxBoundaryCondition` for the gas exchange between carbon dioxide dissolved in the water +specified by the `carbon_chemisty` model, and `air_concentration` with `transfer_velocity` (see +`GasExchangeBoundaryCondition` for details). + +`silicate_and_phosphate_names` should either be `nothing`, a `Tuple`` of symbols specifying the name of the silicate +and phosphate tracers, or a `NamedTuple` of values for the `carbon_chemistry` model. + +`kwargs` are passed on to `GasExchangeBoundaryCondition`. + +Note: The model always requires `T`, `S`, `DIC`, and `Alk` to be present in the model. +""" +function CarbonDioxideGasExchangeBoundaryCondition(; carbon_chemistry = CarbonChemistry(), + transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = CarbonDioxidePolynomialSchmidtNumber()), + air_concentration = 413, # ppmv + wind_speed = 2, + water_concentration = nothing, + silicate_and_phosphate_names = nothing, + kwargs...) + + if isnothing(water_concentration) + water_concentration = CarbonDioxideConcentration(; carbon_chemistry, silicate_and_phosphate_names) + elseif !isnothing(carbon_chemistry) + @warn "Make sure that the `carbon_chemistry` $(carbon_chemistry) is the same as that in `water_concentration` $(water_concentration) (or set it to `nothing`)" + end + + return GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) +end + +""" + OxygenGasExchangeBoundaryCondition(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()), + water_concentration = OxygenConcentration(), + air_concentration = 9352.7, # mmolO₂/m³ + wind_speed = 2, + kwagrs...) + +Returns a `FluxBoundaryCondition` for the gas exchange between oxygen dissolved in the water +specified by the the `OxygenConcentration` in the base model, and `air_concentration` with `transfer_velocity` +(see `GasExchangeBoundaryCondition` for details). + +`kwargs` are passed on to `GasExchangeBoundaryCondition`. +""" +OxygenGasExchangeBoundaryCondition(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()), + water_concentration = OxygenConcentration(), + air_concentration = 9352.7, # mmolO₂/m³ + wind_speed = 2, + kwargs...) = + GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) + +end # module \ No newline at end of file diff --git a/src/Models/GasExchange/carbon_dioxide_concentration.jl b/src/Models/GasExchange/carbon_dioxide_concentration.jl new file mode 100644 index 000000000..0d42f751b --- /dev/null +++ b/src/Models/GasExchange/carbon_dioxide_concentration.jl @@ -0,0 +1,71 @@ +""" + CarbonDioxideConcentration(; carbon_chemistry, + first_virial_coefficient = PolynomialVirialCoefficientForCarbonDioxide(), + cross_viral_coefficient = CrossVirialCoefficientForCarbonDioxide(), + air_pressue = 1 # atm) + +Converts fCO₂ to partial pressure as per Dickson, A.G., Sabine, C.L. and Christian, J.R. (2007), +Guide to Best Practices for Ocean CO 2 Measurements. PICES Special Publication 3, 191 pp. +""" +@kwdef struct CarbonDioxideConcentration{CC<:CarbonChemistry, FV, CV, AP, SP} + carbon_chemistry :: CC + first_virial_coefficient :: FV = PolynomialVirialCoefficientForCarbonDioxide() + cross_virial_coefficient :: CV = CrossVirialCoefficientForCarbonDioxide() + air_pressure :: AP = 1 # atm +silicate_and_phosphate_names :: SP = nothing +end + +summary(::CarbonDioxideConcentration{CC, FV, CV, AP}) where {CC, FV, CV, AP} = + "`CarbonChemistry` derived partial pressure of CO₂ (pCO₂) {$(nameof(CC)), $(nameof(FV)), $(nameof(CV))}" + +@inline function surface_value(cc::CarbonDioxideConcentration, i, j, grid, clock, model_fields) + DIC = @inbounds model_fields.DIC[i, j, grid.Nz] + Alk = @inbounds model_fields.Alk[i, j, grid.Nz] + + T = @inbounds model_fields.T[i, j, grid.Nz] + S = @inbounds model_fields.S[i, j, grid.Nz] + + silicate, phosphate = silicate_and_phosphate(cc.silicate_and_phosphate_names, model_fields) + + fCO₂ = cc.carbon_chemistry(; DIC, Alk, T, S, silicate, phosphate) + + P = surface_value(cc.air_pressure, i, j, grid, clock) * ATM + Tk = T + 273.15 + + B = cc.first_virial_coefficient(Tk) + δ = cc.cross_virial_coefficient(Tk) + + xCO₂ = fCO₂ / P + + # Experimentally this converged xCO₂ to machine precision + for n = 1:3 + φ = P * exp((B + 2 * (1 - xCO₂)^2 * δ) * P / (GAS_CONSTANT * Tk)) + xCO₂ = fCO₂ / φ + end + + return xCO₂ * P # mol/mol to ppm +end + +@inline silicate_and_phosphate(::Nothing, args...) = (0, 0) +@inline silicate_and_phosphate(vals::NamedTuple, args...) = values(vals) +@inline silicate_and_phosphate(val::Tuple, model_fields) = @inbounds getproperty(model_fields, val[1]), getproperty(model_fields, val[2]) + + +# default values from Dickson et al., 2007 +@kwdef struct PolynomialVirialCoefficientForCarbonDioxide{FT} + a :: FT = -1636.75 + b :: FT = 12.0408 + c :: FT = -3.27957 * 10^-2 + d :: FT = 3.16528 * 10^-5 +end + +@inline (fvc::PolynomialVirialCoefficientForCarbonDioxide)(Tk) = + (fvc.a + fvc.b * Tk + fvc.c * Tk^2 + fvc.d * Tk^3) * 10^-6 # cm² / mol to m² / mol + +# default values from Dickson et al., 2007 +@kwdef struct CrossVirialCoefficientForCarbonDioxide{FT} + a :: FT = 57.7 + b :: FT = -0.118 +end + +@inline (cvc::CrossVirialCoefficientForCarbonDioxide)(Tk) = (cvc.a + cvc.b * Tk) * 10^-6 # cm² / mol to m² / mol \ No newline at end of file diff --git a/src/Models/GasExchange/gas_exchange.jl b/src/Models/GasExchange/gas_exchange.jl new file mode 100644 index 000000000..9d05dd5c0 --- /dev/null +++ b/src/Models/GasExchange/gas_exchange.jl @@ -0,0 +1,45 @@ +""" + GasExchange + +`GasExchange` returns the air-sea flux of a gas betwen `water_concentration` and +`air_concentration` with a `transfer_velocity` computed from the temperature +(provided later), and the `wind_speed`. + +`transfer_velocity` should behave as a function of wind speed and temperature (i.e. +`k(u, T)`), `water_concentration` a function of `c(x, y, t, T, field_dependencies...)`. + +`water_concentration`, `air_concentration` and `wind_speed` can either be numbers, +functions of the form `(x, y, t)`, functions of the form `(i, j, grid, clock, model_fields)` +if `discrete_form` is set to true, or any kind of `Field`. + +`water_concentration` should usually be a `[Tracer]Concentration` where is the name of the +tracer (you will have to build your own if this is not `OxygenConcentration`), +or a `CarbonDioxideConcentration` which diagnoses the partial pressure of CO₂ in the water. +""" +struct GasExchange{WS, TV, WC, AC} <: Function + wind_speed :: WS + transfer_velocity :: TV + water_concentration :: WC + air_concentration :: AC +end + +@inline function (g::GasExchange)(i, j, grid, clock, model_fields) + T = @inbounds model_fields.T[i, j, grid.Nz] + + u₁₀ = surface_value(g.wind_speed, i, j, grid, clock) + + k = g.transfer_velocity(u₁₀, T) + + air_concentration = surface_value(g.air_concentration, i, j, grid, clock) + + water_concentration = surface_value(g.water_concentration, i, j, grid, clock, model_fields) + + return k * (water_concentration - air_concentration) +end + +Adapt.adapt_structure(to, g::GasExchange) = GasExchange(adapt(to, g.wind_speed), + adapt(to, g.transfer_velocity), + adapt(to, g.water_concentration), + adapt(to, g.air_concentration)) + +summary(::GasExchange{WS, TV, WC}) where {WS, TV, WC} = "Air-sea `GasExchange` model for $(nameof(WC))" \ No newline at end of file diff --git a/src/Models/GasExchange/gas_transfer_velocity.jl b/src/Models/GasExchange/gas_transfer_velocity.jl new file mode 100644 index 000000000..e97db5a97 --- /dev/null +++ b/src/Models/GasExchange/gas_transfer_velocity.jl @@ -0,0 +1,133 @@ +module ScaledGasTransferVelocity + +export SchmidtScaledTransferVelocity + +using Oceananigans.Units, Adapt + +using OceanBioME.Models.GasExchangeModel: PolynomialParameterisation + +import Adapt: adapt_structure + +""" + SchmidtScaledTransferVelocity(; schmidt_number, + base_transfer_velocity = Ho06()) + +Returns a model for gas transfer velocity which depends on the `u₁₀`, the 10m-wind, and +`T`emperature. The model is of the typical form: + + k(u₁₀, T) = k₆₆₀(u₁₀) √(660/Sc(T)) + +The `base_transfer_velocity` (k₆₆₀) is typically an empirically derived gas transfer velocity +normalised by the Scmidt number for CO₂ at 20°C (660), and the `schmidt_number` (Sc) is a parameterisation +of the gas specific Schmidt number. +""" +@kwdef struct SchmidtScaledTransferVelocity{KB, SC} + base_transfer_velocity :: KB = Ho06() + schmidt_number :: SC +end + +(k::SchmidtScaledTransferVelocity)(u₁₀, T) = k.base_transfer_velocity(u₁₀) * (k.schmidt_number(T) / 660)^(-1/2) + +Adapt.adapt_structure(to, k::SchmidtScaledTransferVelocity) = SchmidtScaledTransferVelocity(adapt(to, k.base_transfer_velocity), + adapt(to, k.schmidt_number)) + +summary(::SchmidtScaledTransferVelocity{KB, SC}) where {KB, SC} = "SchmidtScaledTransferVelocity{$(nameof(KB)), $(nameof(SC))}" +show(io::IO, k::SchmidtScaledTransferVelocity{KB, SC}) where {KB, SC} = + println(io, summary(k), "\n", + " k = k₆₆₀(u₁₀) √(660/Sc(T)),\n", + " Sc(T): $(nameof(SC)),\n", + " k₆₆₀(u₁₀) : $(nameof(KB))") + +""" + Wanninkhof99(scale_factor = 0.0283 / hour / 100) + +Cubic k₆₆₀ parameterisation of Wanninkhof & McGillis (1999) suitable for +short term, in situ wind products. +""" +Wanninkhof99(scale_factor = 0.0283 / hour / 100) = PolynomialParameterisation{3}((0, 0, 0, scale_factor)) + +""" + Ho06(scale_factor = 0.266 / hour / 100) + +Quadratic k₆₆₀ parameterisation of Ho et al. (2006) suitable for the QuickSCAT satellite and short-term +steady wind product. +""" +Ho06(scale_factor = 0.266 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + Nightingale00(linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) + +Cubic k₆₆₀ parameterisation of Nightingale et al. (2000) suitable for +short term, in situ wind products (?). +""" +Nightingale00(linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) = + PolynomialParameterisation{2}((0, linear, quadratic)) + +""" + McGillis01(constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) + +Cubic k₆₆₀ parameterisation of McGillis et al. (2001) suitable for +short term, in situ wind products. +""" +McGillis01(constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) = + PolynomialParameterisation{3}((constant, 0, 0, cubic)) + +""" + Sweeny07(scale_factor = 0.27 / hour / 100) + +Quadratic k₆₆₀ parameterisation of Sweeny et al. (2007) suitable for the +NCEP/NCAR reanalysis 1 product +""" +Sweeny07(scale_factor = 0.27 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + Wanninkhof09(constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) + +Cubic k₆₆₀ parameterisation of Wanninkhof et al (2009) suitable for the +Cross-Calibrated Multi-Platform (CCMP) Winds product +""" +Wanninkhof09(constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) = + PolynomialParameterisation{3}((constant, linear, quadratic, cubic)) + + +""" + Wanninkhof14(scale_factor = 0.251 / hour / 100) + +Quadratic k₆₆₀ parameterisation of Wanninkhof et al (2014) suitable for the +Cross-Calibrated Multi-Platform (CCMP) Winds product +""" +Wanninkhof14(scale_factor = 0.251 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + ERA5(scale_factor = 0.270875 / hour / 100) + +Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) +for the ERA5 wind product by SeaFlux/Luke Gregor et al. (2023). +""" +ERA5(scale_factor = 0.270875 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + JRA55(scale_factor = 0.2601975 / hour / 100) + +Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) +for the JRA55 wind product by SeaFlux/Luke Gregor et al. (2023). +""" +JRA55(scale_factor = 0.2601975 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + NCEP1(scale_factor = 0.2866424 / hour / 100) + +Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) +for the NCEP1 wind product by SeaFlux/Luke Gregor et al. (2023). +""" +NCEP1(scale_factor = 0.2866424 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +""" + CCMP2(scale_factor = 0.256789 / hour / 100) + +Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) +for the CCMP2 wind product by SeaFlux/Luke Gregor et al. (2023). +""" +CCMP2(scale_factor = 0.256789 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) + +end # module \ No newline at end of file diff --git a/src/Models/GasExchange/generic_parameterisations.jl b/src/Models/GasExchange/generic_parameterisations.jl new file mode 100644 index 000000000..7d4c66f0e --- /dev/null +++ b/src/Models/GasExchange/generic_parameterisations.jl @@ -0,0 +1,39 @@ +struct PolynomialParameterisation{N, C} + coefficients :: C +end + +function PolynomialParameterisation{N}(coefficients) where N + length(coefficients) == N + 1 || + throw(ArgumentError("You must provide N+1 coefficients for an order N polynomial")) + + return PolynomialParameterisation{N, typeof(coefficients)}(coefficients) +end + +# generic implementation +@inline function (p::PolynomialParameterisation{N})(x) where N + y = 0 + for n in 0:N + y += @inbounds p.coefficients[n+1] * x^N + end + return y +end + +# fast implementations for low order +@inline (p::PolynomialParameterisation{0})(x) = @inbounds p.coefficients[1] +@inline (p::PolynomialParameterisation{1})(x) = @inbounds p.coefficients[1] + p.coefficients[2] * x +@inline (p::PolynomialParameterisation{2})(x) = @inbounds p.coefficients[1] + p.coefficients[2] * x + p.coefficients[3] * x^2 + +@inline (p::PolynomialParameterisation{3})(x) = @inbounds (p.coefficients[1] + p.coefficients[2] * x + p.coefficients[3] * x^2 + + p.coefficients[4] * x^3) + +@inline (p::PolynomialParameterisation{4})(x) = @inbounds (p.coefficients[1] + p.coefficients[2] * x + p.coefficients[3] * x^2 + + p.coefficients[4] * x^3 + p.coefficients[5] * x^4) + +@inline (p::PolynomialParameterisation{5})(x) = @inbounds (p.coefficients[1] + p.coefficients[2] * x + p.coefficients[3] * x^2 + + p.coefficients[4] * x^3 + p.coefficients[5] * x^4 + p.coefficients[6] * x^5) + + +summary(::PolynomialParameterisation{N}) where N = "Order $N polynomial parameterisation" +show(io::IO, p::PolynomialParameterisation{N}) where N = + println(io, summary(p), "\n", + " p(x) = Σ{n ∈ Z : [0, $N]}(cₙ xⁿ⁻¹) where c = $(p.coefficients)") \ No newline at end of file diff --git a/src/Models/GasExchange/schmidt_number.jl b/src/Models/GasExchange/schmidt_number.jl new file mode 100644 index 000000000..6cf14fcdd --- /dev/null +++ b/src/Models/GasExchange/schmidt_number.jl @@ -0,0 +1,16 @@ + +""" + CarbonDioxidePolynomialSchmidtNumber(; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) + +Schmidt number parameterisation Wanninkhof, 2014 for sea water +""" +CarbonDioxidePolynomialSchmidtNumber(; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) = + PolynomialParameterisation{4}((a, b, c, d, e)) + +""" + OxygenPolynomialSchmidtNumber(; a = 1953.4, b = - 128.0, c = 3.9918, d = -0.050091) + +Schmidt number parameterisation Wanninkhof, 2014 for sea water +""" +OxygenPolynomialSchmidtNumber(; a = 1920.4, b = -135.6, c = 5.2122, d = -0.10939, e = 0.00093777) = + PolynomialParameterisation{4}((a, b, c, d, e)) \ No newline at end of file diff --git a/src/Models/GasExchange/surface_values.jl b/src/Models/GasExchange/surface_values.jl new file mode 100644 index 000000000..6b8e4b5ed --- /dev/null +++ b/src/Models/GasExchange/surface_values.jl @@ -0,0 +1,32 @@ +using Oceananigans.Fields: AbstractField + +# fallback +@inline surface_value(f, args...) = f + +struct ContinuousSurfaceFunction{F} + func :: F +end + +@inline function surface_value(f::ContinuousSurfaceFunction, i, j, grid, clock, args...) + t = clock.time + + x = xnode(i, j, grid.Nz, grid, Center(), Center(), Center()) + y = ynode(i, j, grid.Nz, grid, Center(), Center(), Center()) + + return f.func(x, y, t, args...) +end + +struct DiscreteSurfaceFuncton{F} + func :: F +end + +@inline surface_value(f::DiscreteSurfaceFuncton, i, j, grid, clock, args...) = f.func(i, j, grid, clock, args...) + +# GPU compatible, mainly for fields +@inline surface_value(f::AbstractArray{<:Any, 2}, i, j, grid, clock, args...) = @inbounds f[i, j] +@inline surface_value(f::AbstractArray{<:Any, 3}, i, j, grid, clock, args...) = @inbounds f[i, j, grid.Nz] + +# fallback +normalise_surface_function(f; kwargs...) = f + +normalise_surface_function(f::Function; discrete_form = false) = discrete_form ? DiscreteSurfaceFuncton(f) : ContinuousSurfaceFunction(f) \ No newline at end of file diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 6d8a719b4..72f389ab1 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -18,6 +18,8 @@ Optional """ module SLatissimaModel +export SLatissima + using Roots, KernelAbstractions using OceanBioME.Particles: BiogeochemicalParticles, get_node using Oceananigans.Units diff --git a/src/Models/Models.jl b/src/Models/Models.jl index c50245e38..f7269e855 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -1,3 +1,37 @@ +module Models + +export Sediments + +export NPZD, + NutrientPhytoplanktonZooplanktonDetritus, + LOBSTER + +export SLatissima + +export CarbonChemistry + +export GasExchange, + CarbonDioxideGasExchangeBoundaryCondition, + OxygenGasExchangeBoundaryCondition, + GasExchangeBoundaryCondition, + ScaledGasTransferVelocity, + SchmidtScaledTransferVelocity, + CarbonDioxidePolynomialSchmidtNumber, + OxygenPolynomialSchmidtNumber + +include("Sediments/Sediments.jl") include("AdvectedPopulations/LOBSTER/LOBSTER.jl") include("AdvectedPopulations/NPZD.jl") -include("Individuals/SLatissima.jl") \ No newline at end of file +include("Individuals/SLatissima.jl") +include("seawater_density.jl") +include("CarbonChemistry/CarbonChemistry.jl") +include("GasExchange/GasExchange.jl") + +using .Sediments +using .LOBSTERModel +using .NPZDModel +using .SLatissimaModel +using .CarbonChemistryModel +using .GasExchangeModel + +end # module \ No newline at end of file diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Models/Sediments/Sediments.jl similarity index 100% rename from src/Boundaries/Sediments/Sediments.jl rename to src/Models/Sediments/Sediments.jl diff --git a/src/Boundaries/Sediments/coupled_timesteppers.jl b/src/Models/Sediments/coupled_timesteppers.jl similarity index 98% rename from src/Boundaries/Sediments/coupled_timesteppers.jl rename to src/Models/Sediments/coupled_timesteppers.jl index 63c3ffd87..a603dede4 100644 --- a/src/Boundaries/Sediments/coupled_timesteppers.jl +++ b/src/Models/Sediments/coupled_timesteppers.jl @@ -1,5 +1,5 @@ using Oceananigans: NonhydrostaticModel, prognostic_fields, HydrostaticFreeSurfaceModel -using OceanBioME.Boundaries.Sediments: AbstractSediment +using OceanBioME.Models.Sediments: AbstractSediment using Oceananigans.TimeSteppers: ab2_step_field!, rk3_substep_field!, stage_Δt using Oceananigans.Utils: work_layout, launch! using Oceananigans.TurbulenceClosures: implicit_step! diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Models/Sediments/instant_remineralization.jl similarity index 100% rename from src/Boundaries/Sediments/instant_remineralization.jl rename to src/Models/Sediments/instant_remineralization.jl diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Models/Sediments/simple_multi_G.jl similarity index 98% rename from src/Boundaries/Sediments/simple_multi_G.jl rename to src/Models/Sediments/simple_multi_G.jl index b9b64ac97..138049c27 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Models/Sediments/simple_multi_G.jl @@ -76,14 +76,14 @@ so it will therefore be important to also thoroughly validate the oxygen model ( Example ======= -```jldoctest simplemultig; filter = r".*@ OceanBioME.Boundaries.Sediments.*" +```jldoctest simplemultig; filter = r".*@ OceanBioME.Models.Sediments.*" julia> using OceanBioME, Oceananigans, OceanBioME.Sediments julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200)); julia> sediment_model = SimpleMultiG(; grid) ┌ Warning: Sediment models are an experimental feature and have not yet been validated. -└ @ OceanBioME.Boundaries.Sediments ~/OceanBioME.jl/src/Boundaries/Sediments/simple_multi_G.jl:102 +└ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:104 [ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC. Single-layer multi-G sediment model (Float64) ``` diff --git a/src/Models/seawater_density.jl b/src/Models/seawater_density.jl new file mode 100644 index 000000000..650638187 --- /dev/null +++ b/src/Models/seawater_density.jl @@ -0,0 +1,39 @@ +using GibbsSeaWater: gsw_rho, gsw_sa_from_sp +using SeawaterPolynomials.TEOS10: _ρ, τ, s, ζ + +""" + teos10_density(T, Sp, Pbar = 0, lon = 0, lat = 0) + +Returns sea water density computed by `GibbsSeaWater.jl` from +`T` in degrees centigrade, `Sp` in PSU, `P`ressure in bar. Location +`lat`itude and `lon`gitude are required for accurate conversion from +practical to absolute salinity. + +This function *will not* run on GPU so can not be used as the default. +""" +@inline function teos10_density(T, Sp, Pbar = 0, lon = 0, lat = 0) + Pdbar = 10 * Pbar + + Sa = gsw_sa_from_sp(Sp, Pdbar, lon, lat) + + return gsw_rho(Sa, T, Pdbar) +end + +""" + teos10_polynomial_approximation(T, Sp, Pbar = 0, lon = 0, lat = 0) + +Returns sea water density computed by `SeawaterPolynomials.jl` TEOS10 +polynomial approximation from `T` in degrees centigrade, `S` in PSU, +`P`ressure in bar. This approximation also requires us to approximate +geopotential depth from pressure at 1 m / dbar, and absolute salinity +as practial salinity. + +This method is less accurate than `teos10_density` but will run on GPU. +""" +@inline function teos10_polynomial_approximation(T, Sp, Pbar = 0, lon = 0, lat = 0) + Z = 10 * Pbar # 10 m / bar + + Sa = Sp # error is *very* small in the pH and pCO₂ from this approximation + + return _ρ(τ(T), s(Sa), ζ(Z)) +end \ No newline at end of file diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index fabea3538..3f8765ca6 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -19,8 +19,14 @@ export Particles # Light models export TwoBandPhotosyntheticallyActiveRadiation, PrescribedPhotosyntheticallyActiveRadiation -# Boundaries -export Boundaries, Sediments, GasExchange, FlatSediment +# airsea flux +export GasExchange, CarbonDioxideGasExchangeBoundaryCondition, OxygenGasExchangeBoundaryCondition, GasExchangeBoundaryCondition + +# carbon chemistry +export CarbonChemistry + +# sediment +export Sediments, FlatSediment # Utilities export column_advection_timescale, sinking_advection_timescale, Budget @@ -175,17 +181,13 @@ show(io::IO, model::Biogeochemistry) = summary(modifiers::Tuple) = tuple([summary(modifier) for modifier in modifiers]) include("Utils/Utils.jl") -include("Boundaries/Boundaries.jl") include("Light/Light.jl") include("Particles/Particles.jl") include("BoxModel/boxmodel.jl") include("Models/Models.jl") -using .Boundaries using .Light using .BoxModels -using .LOBSTERModel -using .NPZDModel -import .SLatissimaModel.SLatissima +using .Models end #module diff --git a/test/runtests.jl b/test/runtests.jl index 04cad701e..44cde13f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ include("test_light.jl") include("test_slatissima.jl") include("test_LOBSTER.jl") include("test_NPZD.jl") -include("test_gasexchange.jl") +include("test_gasexchange_carbon_chem.jl") include("test_sediments.jl") if architecture == CPU() diff --git a/test/test_gasexchange.jl b/test/test_gasexchange.jl deleted file mode 100644 index eb3347de0..000000000 --- a/test/test_gasexchange.jl +++ /dev/null @@ -1,69 +0,0 @@ -include("dependencies_for_runtests.jl") - -using OceanBioME: Boundaries, GasExchange, LOBSTER -using Oceananigans, DataDeps, JLD2, Statistics -using Oceananigans.Units - -const year = years = 365days # just for the idealised case below - -dd = DataDep( - "test_data", - "CODAP-NA (https://essd.copernicus.org/articles/13/2777/2021/) data for testing pCO₂ calculations", - "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/CODAP_data.jld2" -) - -register(dd) - -function test_gas_exchange_model(grid, air_concentration) - model = NonhydrostaticModel(; grid, - tracers = (:T, :S), - biogeochemistry = LOBSTER(; grid, carbonates = true), - boundary_conditions = (DIC = FieldBoundaryConditions(top = GasExchange(; air_concentration, gas = :CO₂)), )) - - @test isa(model.tracers.DIC.boundary_conditions.top.condition.func, GasExchange) - - set!(model, T = 15.0, S = 35.0, DIC = 2220, Alk = 2500) - - # is everything communicating properly? (can't think of a way to not use allow scalar here) - value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) - - @test value ≈ -0.0003 atol = 0.0001 - - # just incase we broke something - @test isnothing(time_step!(model, 1.0)) - - return nothing -end - -@testset "Gas exchange values" begin - # approximatly correct sized pCO₂ from DIC and ALK - pCO₂_model = Boundaries.OCMIP_default - - @load datadep"test_data/CODAP_data.jld2" DIC Alk T S pH pCO₂ - - pCO₂_results = similar(pCO₂) - - for (idx, DIC) in enumerate(DIC) - pCO₂_results[idx] = pCO₂_model(DIC, Alk[idx], T[idx], S[idx]) - end - - pCO₂_err = pCO₂ .- pCO₂_results - - # not great - @test (mean(pCO₂_err) < 10 && std(pCO₂_err) < 15) -end - -grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) - -@inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) - -#conc_field = CenterField(grid, indices=(:, :, grid.Nz)) -#set!(conc_field, 413.0) -#Oceananigans.fill_halo_regions!(conc_field) - -@testset "Gas exchange coupling" begin - for air_concentration in [413.1, conc_function]#, conc_field] - @info "Testing with $(typeof(air_concentration))" - test_gas_exchange_model(grid, air_concentration) - end -end \ No newline at end of file diff --git a/test/test_gasexchange_carbon_chem.jl b/test/test_gasexchange_carbon_chem.jl new file mode 100644 index 000000000..b1642ad55 --- /dev/null +++ b/test/test_gasexchange_carbon_chem.jl @@ -0,0 +1,137 @@ +include("dependencies_for_runtests.jl") + +using Oceananigans, Oceananigans.Units, DataDeps, JLD2, Statistics + +using Oceananigans.Fields: ConstantField + +using OceanBioME: GasExchange, LOBSTER, CarbonChemistry +using OceanBioME.Models.GasExchangeModel: surface_value + +using OceanBioME.Models.CarbonChemistryModel: IonicStrength, K0, K1, K2, KB, KW, KS, KF, KP, KSi, KSP_aragonite, KSP_calcite + +const year = years = 365days # just for the idealised case below + +dd = DataDep( + "test_data", + "CODAP-NA (https://essd.copernicus.org/articles/13/2777/2021/) data for testing pCO₂ calculations", + "https://github.com/OceanBioME/OceanBioME_example_data/raw/main/CODAP_data.jld2" +) + +register(dd) + +function test_gas_exchange_model(grid, air_concentration) + model = NonhydrostaticModel(; grid, + tracers = (:T, :S), + biogeochemistry = LOBSTER(; grid, carbonates = true), + boundary_conditions = (DIC = FieldBoundaryConditions(top = CarbonDioxideGasExchangeBoundaryCondition(; air_concentration)), )) + + set!(model, T = 15.0, S = 35.0, DIC = 2220, Alk = 2500) + + # is everything communicating properly? (can't think of a way to not use allow scalar here) + value = CUDA.@allowscalar Oceananigans.getbc(model.tracers.DIC.boundary_conditions.top, 1, 1, grid, model.clock, fields(model)) + + return isa(model.tracers.DIC.boundary_conditions.top.condition.func, GasExchange)&&≈(value, -0.0002; atol = 0.0001)&&isnothing(time_step!(model, 1.0)) +end + +@testset "pCO₂ values" begin + # approximatly correct sized pCO₂ from DIC and Alk + fCO₂_model = CarbonChemistry() + + @load datadep"test_data/CODAP_data.jld2" DIC Alk T S pH pCO₂ + + fCO₂_results = similar(pCO₂) + pH_results = similar(pH) + + for (idx, DIC) in enumerate(DIC) + fCO₂_results[idx] = fCO₂_model(; DIC, Alk = Alk[idx], T = T[idx], S = S[idx]) + pH_results[idx] = fCO₂_model(; DIC, Alk = Alk[idx], T = T[idx], S = S[idx], return_pH = true) + end + + fCO₂_err = pCO₂ .- fCO₂_results + pH_err = pH .- pH_results + + # not great not terrible + @test (mean(fCO₂_err) < 9 && std(fCO₂_err) < 10) && (mean(pH_err) < 0.01 && std(pH_err) < 0.01) +end + +grid = RectilinearGrid(architecture; size=(1, 1, 2), extent=(1, 1, 1)) + +@inline conc_function(x, y, t) = 413.0 + 10.0 * sin(t * π / year) + +conc_field = CenterField(grid) + +set!(conc_field, (args...) -> 10 * randn() + 413) + +@testset "Gas exchange coupling" begin + for air_concentration in [413.1, conc_function, conc_field] + @info "Testing gas exchange with $(typeof(air_concentration))" + @test test_gas_exchange_model(grid, air_concentration) + end +end + +@testset "Carbon chemistry" begin + carbon_chemistry = CarbonChemistry() + + # test conditions + S = 35 + Tk = 298.15 + P = 300 + + # values from Dickson et. al, 2007 + @test ≈(log(carbon_chemistry.solubility(Tk, S)), -3.5617; atol=0.0001) + @test ≈(log10(carbon_chemistry.carbonic_acid.K1(Tk, S)), -5.8472; atol=0.0001) + @test ≈(log10(carbon_chemistry.carbonic_acid.K2(Tk, S)), -8.9660; atol=0.0001) + @test ≈(log(carbon_chemistry.boric_acid(Tk, S)), -19.7964; atol=0.0001) + @test ≈(log(carbon_chemistry.water(Tk, S)), -30.434; atol=0.001) + @test ≈(log(carbon_chemistry.sulfate(Tk, S)), -2.30; atol=0.01) + @test ≈(log(carbon_chemistry.fluoride(Tk, S)), -6.09; atol=0.01) + @test ≈(log(carbon_chemistry.phosphoric_acid.KP1(Tk, S)), -3.71; atol=0.01) + @test ≈(log(carbon_chemistry.phosphoric_acid.KP2(Tk, S)), -13.727; atol=0.001) + @test ≈(log(carbon_chemistry.phosphoric_acid.KP3(Tk, S)), -20.24; atol=0.01) + @test ≈(log(carbon_chemistry.silicic_acid(Tk, S)), -21.61; atol=0.01) + + # values from Zeebe & Wolf-Gladrow, 2001 + @test ≈(carbon_chemistry.carbonic_acid.K1.pressure_correction(Tk, P), 1.30804; atol=0.00001) + @test ≈(carbon_chemistry.carbonic_acid.K2.pressure_correction(Tk, P), 1.21341; atol=0.00001) + @test ≈(carbon_chemistry.boric_acid.pressure_correction(Tk, P), 1.38024; atol=0.00001) + @test ≈(carbon_chemistry.water.pressure_correction(Tk, P), 1.23784; atol=0.00001) + @test ≈(carbon_chemistry.sulfate.pressure_correction(Tk, P), 1.21844; atol=0.00001) + @test ≈(carbon_chemistry.fluoride.pressure_correction(Tk, P), 1.13151; atol=0.00001) + @test ≈(carbon_chemistry.phosphoric_acid.KP1.pressure_correction(Tk, P), 1.14852; atol=0.00001) + @test ≈(carbon_chemistry.phosphoric_acid.KP2.pressure_correction(Tk, P), 1.27298; atol=0.00001) + @test ≈(carbon_chemistry.phosphoric_acid.KP3.pressure_correction(Tk, P), 1.32217; atol=0.00001) + + # Calcite and aragonite solubility + KspA = KSP_aragonite() + KspC = KSP_calcite() + + # Zeebe & Wolf-Gladrow, 2001, Appendix A + @test ≈(log(KspA(Tk, S)), -6.1883; atol = 0.0001) + @test ≈(log(KspC(Tk, S)), -6.3693; atol = 0.0001) + + @test ≈(KspA.pressure_correction(Tk, P), 1.47866; atol=0.00001) + @test ≈(KspC.pressure_correction(Tk, P), 1.52962; atol=0.00001) +end + +@testset "Gas exchange constants defaults" begin + CO₂_exchange = CarbonDioxideGasExchangeBoundaryCondition().condition.func + O₂_exchange = OxygenGasExchangeBoundaryCondition().condition.func + + T = 20 + + # values from Wanninkhof, 2014 + @test ≈(CO₂_exchange.transfer_velocity.schmidt_number(T), 668, atol = 1) + @test ≈(O₂_exchange.transfer_velocity.schmidt_number(T), 568, atol = 1) + + Tk = 25+273.15 + # values from Dickson et. al, 2007 + @test ≈(CO₂_exchange.water_concentration.first_virial_coefficient(Tk), -123.2 * 10^-6, atol=10^-8) + @test ≈(CO₂_exchange.water_concentration.cross_virial_coefficient(Tk), 22.5 * 10^-6, atol=10^-7) + + T = ConstantField(25) + S = ConstantField(35) + DIC = ConstantField(2136.242890518708) + Alk = ConstantField(2500) + # value from Dickson et. al, 2007 + @test ≈(surface_value(CO₂_exchange.water_concentration, 1, 1, BoxModelGrid(), Clock(; time = 0), (; T, S, DIC, Alk)), 350, atol = 0.1) +end \ No newline at end of file diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 503ef93f5..cda5478f7 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -9,7 +9,7 @@ using Oceananigans.Fields: TracerFields using Oceananigans.Operators: volume, Azᶠᶜᶜ -using OceanBioME.LOBSTERModel: VariableRedfieldLobster +using OceanBioME.Models.LOBSTERModel: VariableRedfieldLobster function intercept_tracer_tendencies!(model, intercepted_tendencies) for (name, field) in enumerate(intercepted_tendencies) diff --git a/test/test_slatissima.jl b/test/test_slatissima.jl index a9aed5c0e..071fd47ee 100644 --- a/test/test_slatissima.jl +++ b/test/test_slatissima.jl @@ -1,6 +1,6 @@ include("dependencies_for_runtests.jl") -using OceanBioME.SLatissimaModel: SLatissima +using OceanBioME: SLatissima using Oceananigans.Units using Oceananigans.Fields: TracerFields using Oceananigans.Architectures: on_architecture diff --git a/validation/glodap_fit.jl b/validation/glodap_fit.jl new file mode 100644 index 000000000..1eb7e4660 --- /dev/null +++ b/validation/glodap_fit.jl @@ -0,0 +1,144 @@ +using CSV, DataFrames, Random, CairoMakie +using OceanBioME: CarbonChemistry + +using OceanBioME.Models: teos10_density, teos10_polynomial_approximation +using OceanBioME.Models.CarbonChemistryModel: K0, K1, K2, KF + + + +# get the glodap data - downloaded "merged" file from https://www.ncei.noaa.gov/data/oceans/ncei/ocads/data/0283442/ +data = CSV.read("validation/GLODAP.csv", DataFrame) + +DIC_name = "G2tco2" +Alk_name = "G2talk" +T_name = "G2temperature" +S_name = "G2salinity" +P_name = "G2pressure" +Si_name = "G2silicate" +PO_name = "G2phosphate" +pH_name = "G2phtsinsitutp" +depth_name = "G2depth" +cruise_name = "G2cruise" +fCO2_name = "G2fco2" +lon_name = "G2longitude" +lat_name = "G2latitude" + +DIC_name, Alk_name, T_name, S_name, P_name, Si_name, PO_name, pH_name + +# remove low quality data and non QCed +# and drop a weird point from cruise with a massive errors (only 2 points) +#=data = filter(row -> any([row[name*"f"] == 2 for name in [DIC_name, Alk_name, S_name]]) && + !any([row[name] == -9999 for name in [DIC_name, Alk_name, T_name, S_name]]) && + #row[cruise_name] != 5005 && + #row[cruise_name] != 345 && + row[Alk_name] > 1500 && + row["G2phtsqc"] == 1, data) +=# +# if we don't have Si and PO₄ the error gets a lot bigger but not in the surface values (as much) +data = filter(row -> any([row[name*"f"] == 2 for name in [DIC_name, Alk_name, S_name, pH_name]]) && + !any([row[name] == -9999 for name in [DIC_name, Alk_name, T_name, S_name, P_name, Si_name, PO_name, pH_name]]) && + row["G2expocode"] != "325020210316" && + row["G2expocode"] != "33RO20071215" && + row[cruise_name] != 270 && + #row[Alk_name] > 1500 && + row["G2phtsqc"] == 1, + data) + +Nrows = size(data, 1) + +# setup the model + +density_function = teos10_density + +# pH_error = 0.00022 ± 0.01223 -> mean is zero +# pCO₂_error = -3 ± 28.8 -> mean is also zero +# with teos polynomial appriximation error is almost identical in pCO2 but slightly bigger in pH +carbon_chemistry = CarbonChemistry() + +# compare +glodap_pH = zeros(Nrows) +glodap_pCO₂ = zeros(Nrows) + +model_pH = zeros(Nrows) +model_pCO₂ = zeros(Nrows) + +pH_error = zeros(Nrows) +pCO₂_error = zeros(Nrows) + +Threads.@threads for n in 1:Nrows + T = data[n, T_name] + S = data[n, S_name] + P = data[n, P_name] * 0.1 # dbar to bar + + lon = data[n, lon_name] + lat = data[n, lat_name] + + density = density_function(T, S, P, lon, lat) + + DIC = data[n, DIC_name] * density * 1e-3 + Alk = data[n, Alk_name] * density * 1e-3 + + silicate = ifelse(silicate == -9999, 0, silicate) + phosphate = ifelse(phosphate == -9999, 0, phosphate) + + silicate = data[n, Si_name] * density * 1e-3 + phosphate = data[n, PO_name] * density * 1e-3 + + P = ifelse(P == -9999 * 0.1, 0, P) + + if data[n, pH_name*"f"] == 2 && data[n, pH_name] != -9999 + model_pH[n] = carbon_chemistry(; DIC, Alk, T, S, P, silicate, phosphate, lon, lat, return_pH = true) + glodap_pH[n] = data[n, pH_name] + pH_error[n] = model_pH[n] - glodap_pH[n] + else + model_pH[n] = NaN + glodap_pH[n] = NaN + pH_error[n] = NaN + end + + if data[n, fCO2_name*"f"] == 2 && data[n, fCO2_name*"temp"] != -9999 + # fCO2 is reported at no pressure and at a different temp + model_pCO₂[n] = carbon_chemistry(; DIC, Alk, T = data[n, fCO2_name*"temp"], S, silicate, phosphate, lon, lat) + glodap_pCO₂[n] = data[n, fCO2_name] + pCO₂_error[n] = model_pCO₂[n] - glodap_pCO₂[n] + else + model_pCO₂[n] = NaN + glodap_pCO₂[n] = NaN + pCO₂_error[n] = NaN + end + +end + +function plot_errors(; subset = :, + to_plot = pH_error, + title = "ΔpH", + color_name = S_name, + color_marker = data[:, color_name][subset], + limit_y = false, + ylim = 100) + fig = Figure(); + + ax = Axis(fig[1, 1], title = "DIC")#, aspect = DataAspect(), title = "Local grid") + ax2 = Axis(fig[1, 2], title = "Alk") + ax3 = Axis(fig[1, 3], title = "pH") + ax4 = Axis(fig[2, 1], title = "T") + ax5 = Axis(fig[2, 2], title = "S") + ax6 = Axis(fig[2, 3], title = "Depth") + + sc=scatter!(ax, data[:, DIC_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + scatter!(ax2, data[:, Alk_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + scatter!(ax3, data[:, pH_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + scatter!(ax4, data[:, T_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + scatter!(ax5, data[:, S_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + scatter!(ax6, data[:, depth_name][subset], to_plot[subset], markersize = 3, alpha = 0.5, color = color_marker) + + if limit_y + [ylims!(a, -ylim, ylim) for a in (ax, ax2, ax3, ax4, ax5, ax6)] + end + + Colorbar(fig[1:2, 4], sc) + + Label(fig[0, :], title) + + return fig +end \ No newline at end of file