Library

A fast and flexible modelling environment for modelling the coupled interactions between ocean biogeochemistry, carbonate chemistry, and physics.

(day_length::CBMDayLength)(t, φ)

Returns the length of day in seconds at the latitude φ, t seconds after the start of the year.

ScaleNegativeTracers(; tracers, scalefactors = ones(length(tracers)), warn = false, invalid_fill_value = NaN)

Constructs a modifier to scale tracers so that none are negative. Use like:

modifier = ScaleNegativeTracers((:P, :Z, :N))
biogeochemistry = Biogeochemistry(...; modifier)

This method is better, though still imperfect, method to prevent numerical errors that lead to negative tracer values compared to ZeroNegativeTracers. Please see discussion in github.

Future plans include implement a positivity-preserving timestepping scheme as the ideal alternative.

~~If warn is true then scaling will raise a warning.~~

invalid_fill_value specifies the value to set the total cell content to if the total is less than 0 (meaning that total tracer conservation can not be enforced). If the value is set to anything other than NaN this scheme no longer conserves mass. While this may be useful to prevent spurious numerics leading to crashing care should be taken that the mass doesn't deviate too much.

This scheme is similar to that used by NEMO-PISCES, although they scale the tendency rather than the value, while other Earth system models simply set negative tracers to zero, for example NCAR's MARBL and NEMO-TOPAZ2, which does not conserve mass. More complicated schemes exist, for example ROMS-BECS uses an implicite-itterative approach where each component is updated in sequence to garantee mass conservation, possibly at the expense of numerical precision.

ScaleNegativeTracers(bgc::AbstractBiogeochemistry; warn = false)

Construct a modifier to scale the conserved tracers in bgc biogeochemistry.

If warn is true then scaling will raise a warning.

ZeroNegativeTracers(; exclude = ())

Construct a modifier that zeroes any negative tracers excluding those listed in exclude.

Biogeochemistry(underlying_biogeochemistry;
                light_attenuation = nothing,
                sediment = nothing,
                particles = nothing,
                modifiers = nothing)

Construct a biogeochemical model based on underlying_biogeochemistry which may have a light_attenuation model, sediment, particles, and modifiers.

Keyword Arguments

• light_attenuation: light attenuation model which integrated the attenuation of available light
• sediment_model: slot for AbstractSediment
• particles: slot for BiogeochemicalParticles
• modifiers: slot for components which modify the biogeochemistry when the tendencies have been calculated or when the state is updated

conserved_tracers(model::UnderlyingBiogeochemicalModel, args...; kwargs...)

Returns the names of tracers which together are conserved in model

redfield(i, j, k, val_tracer_name, bgc, tracers)

Returns the redfield ratio of tracer_name from bgc at i, j, k.

redfield(val_tracer_name, bgc)
redfield(val_tracer_name, bgc, tracers)

Returns the redfield ratio of tracer_name from bgc when it is constant across the domain.

Nutrient-Phytoplankton-Zooplankton-Detritus model of Kuhn et al. (2015).

Tracers

• Nutrients: N (mmol N/m³)
• Phytoplankton: P (mmol N/m³)
• Zooplankton: Z (mmol N/m³)
• Detritus: D (mmol N/m³)

Required submodels

• Photosynthetically available radiation: PAR (W/m²)

NutrientPhytoplanktonZooplanktonDetritus(; grid::AbstractGrid{FT},
                                           initial_photosynthetic_slope::FT = 0.1953 / day, # 1/(W/m²)/s
                                           base_maximum_growth::FT = 0.6989 / day, # 1/s
                                           nutrient_half_saturation::FT = 2.3868, # mmol N/m³
                                           base_respiration_rate::FT = 0.066 / day, # 1/s/(mmol N / m³)
                                           phyto_base_mortality_rate::FT = 0.0101 / day, # 1/s/(mmol N / m³)
                                           maximum_grazing_rate::FT = 2.1522 / day, # 1/s
                                           grazing_half_saturation::FT = 0.5573, # mmol N/m³
                                           assimulation_efficiency::FT = 0.9116, 
                                           base_excretion_rate::FT = 0.0102 / day, # 1/s/(mmol N / m³)
                                           zoo_base_mortality_rate::FT = 0.3395 / day, # 1/s/(mmol N / m³)²
                                           remineralization_rate::FT = 0.1213 / day, # 1/s

                                           surface_photosynthetically_active_radiation = default_surface_PAR,
                                           light_attenuation::LA =
                                               TwoBandPhotosyntheticallyActiveRadiation(; grid,
                                                                                          surface_PAR = surface_photosynthetically_active_radiation),
                                           sediment::S = nothing,
            
                                           sinking_speeds = (P = 0.2551/day, D = 2.7489/day),
                                           open_bottom::Bool = true,

                                           scale_negatives = false,
                                                                                  
                                           particles::P = nothing,
                                           modifiers::M = nothing)

Construct a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) biogeochemical model.

Keyword Arguments

• grid: (required) the geometry to build the model on, required to calculate sinking
• initial_photosynthetic_slope, ..., remineralization_rate: NPZD parameter values
• surface_photosynthetically_active_radiation: function (or array in the future) for the photosynthetically available radiation at the surface, should be shape f(x, y, t)
• light_attenuation: light attenuation model which integrated the attenuation of available light
• sediment: slot for BiogeochemicalSediment
• sinking_speed: named tuple of constant sinking, of fields (i.e. ZFaceField(...)) for any tracers which sink (convention is that a sinking speed is positive, but a field will need to follow the usual down being negative)
• open_bottom: should the sinking velocity be smoothly brought to zero at the bottom to prevent the tracers leaving the domain
• scale_negatives: scale negative tracers?
• particles: slot for BiogeochemicalParticles
• modifiers: slot for components which modify the biogeochemistry when the tendencies have been calculated or when the state is updated

Example

julia> using OceanBioME

julia> using Oceananigans

julia> grid = RectilinearGrid(size=(20, 30), extent=(200, 200), topology=(Bounded, Flat, Bounded));

julia> model = NutrientPhytoplanktonZooplanktonDetritus(; grid)
NutrientPhytoplanktonZooplanktonDetritus{Float64} model, with (:P, :D) sinking 
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: Nothing
 Particles: Nothing
 Modifiers: Nothing

The Lodyc-DAMTP Ocean Biogeochemical Simulation Tools for Ecosystem and Resources (LOBSTER) model.

LOBSTER is an NPZD type model with multiple nutrients and detritus pools. In the default configuration there is one phytoplankton and one zooplankton, nitrate and ammonia, and small, large, and dissovled detritus.

Optionally it can also include iron (which limits phytoplankton growth), the carbonate system (DIC and alkalinity), oxygen, and variable N:C ratio detritus where each of dissolved, small, and large particles have both nitrogen and carbon compartements.

Required submodels

• Photosynthetically available radiation: PAR (W/m²)

When carbonate system is active:

• Temperature: T (ᵒC)
• Salinity: S (‰)

CarbonateSystem

CarbonateSystem defines the carbonate system for the LOBSTER biogeochemical model and evolves dissolved inorganic carbon (DIC) and alkalinity (Alk). DIC is the total inorganic carbon, typically considered to be dissolved carbon dioxide, carbonic acid, bicarbonate, and carbonate. Alkalinity is the buffering capacity of the water which is approximatly the amount of bases in solution that can be neutralised by the addition of acid so is approximatly the sum of non-conserved ions such as carbonate, bicarbonate, hydroxide ions, etc.

With this module activated the carbon budget can be closed as the flux of carbon into the biology and detritus comes from the DIC. Additionally, the flux in/out of the water from the air can be computed via the CarbonChemistry module.

DIC is produced by remineralisation of organic matter, for example the breakdown of detritus or excretion from zooplankton waste. It is also produced by the dissolution of calcite in the zooplankton gut. It is removed by photosynthesis in phytoplankton.

Alkalinity is produced (by the removal of nitrose acid (?)), and removed by the uptake of calcite into phytoplankton.

DIC and Alk concentration is only one way coupled with the rest of the biogeochemistry and does not effect any other groups (e.g. acidifcation does not effect phytoplankton growth). To capture this effect a different biology could be defined.

LOBSTER(; grid

          nutrients = NitrateAmmonia(),
          biology = PhytoZoo(),
          detritus = TwoParticleAndDissolved(grid),
          carbonate_system = nothing,
          oxygen = nothing,

          surface_photosynthetically_active_radiation = default_surface_PAR,

          light_attenuation =
                TwoBandPhotosyntheticallyActiveRadiation(; grid, 
                                                           surface_PAR = surface_photosynthetically_active_radiation),
            
          sediment = nothing,

          scale_negatives = false,
          invalid_fill_value = NaN,

          particles = nothing,
          modifiers = nothing)

Construct an instance of the LOBSTER biogeochemical model.

Keyword Arguments

• grid: (required) the geometry to build the model on, required to configure sinking speeds
• nutrients: nutrint component, defaults to NitrateAmmonia which has nitrate (NO₃) and ammonia (NH₄),

currently NitrateAmmoniaIron is the alternative which includes iron (Fe)

• biology: biological (living) component, defaults to PhytoZoo which has phytoplankton (P) and zooplankton (Z)
• detritus: non-living organic component, default to TwoParticlesAndDissolved which has small and large

particles (sPOM and bPOM) and dissolved organic (DOM), alternativly can be VariableRedfieldDetritus which has sPON, bPON, DON, sPOC, bPOC, and DOC

• carbonate_system: an optional model component for the inorganic carbon components, defaults to nothing but

can be CarbonateSystem which has dissolved inorganic carbon (DIC) and alkalinity (Alk)

• oxygen: an optional model component for oxygen, defaults to nothing but can be Oxygen which has oxygen (O₂)
• surface_photosynthetically_active_radiation: funciton for the photosynthetically available radiation at the surface, should be shape f(x, y, t)
• light_attenuation: light attenuation model which integrated the attenuation of available light
• sediment_model: slot for BiogeochemicalSediment
• scale_negatives: scale negative tracers?
• particles: slot for BiogeochemicalParticles
• modifiers: slot for components which modify the biogeochemistry when the tendencies have been calculated or when the state is updated

Example

julia> using OceanBioME, Oceananigans

julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));

julia> model = LOBSTER(; grid)
LOBSTER model (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM) 
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: Nothing
 Particles: Nothing
 Modifiers: Nothing

NitrateAmmonia

NitrateAmmonia defines the default nutrient component for the LOBSTER biogeochemical model which includes nitrate (NO₃), and ammonia (NH₄) in mmol N / m³.

Nitrate is only taken up by the biological component (by default only the phytoplankton), and replenished by the nitrifcation of ammonia at the nitrification_rate.

Ammonia is also taken up by the biological component, and lost as it is turned into ammonia. It is replaced through waste from the biological system (phytoplankton exudation, and zooplankton excretion), and from the degredation of detritus.

NitrateAmmoniaIron

NitrateAmmoniaIron defines a nutrient component for the LOBSTER biogeochemical model which includes nitrate (NO₃), and ammonia (NH₄) in mmol N / m³ in the same way as the NitrateAmmonia component, but with the addition of an iron component (Fe) in mmol Fe / m³.

Iron is only taken up by biological component and not replenished, as iron gets remineralised a lot slower than the nitrogen and carbon in detritus so the iron is "never" returned to the biologically available pool. Instead it is usually replenished by surface deposition of dust. For example, Signorini, et al. 2001 uses the same formulation.

Signorini, S. R., C. R. McClain, J. R. Christian, and C. S. Wong (2001), Seasonal and interannual variability of phytoplankton, nutrients, TCO2, pCO2, and O2 in the eastern subarctic Pacific (ocean weather station Papa), J. Geophys. Res., 106(C12), 31197–31215, doi:10.1029/2000JC000343.

Oxygen

Oxygen defines the evolution of oxygen (O₂) concentration for the LOBSTER biogeochemical model.

Oxygen is produced by photosynthesis in phytoplankton, and removed by nitrate production and oxidation of ammonia.

Oxygen concentration is only one way coupled with the rest of the biogeochemistry and does not effect any other groups (e.g. low oxygen does not reduce zooplankton growth). To capture this effect a different biology could be defined.

PhytoZoo

PhytoZoo defines the default living component for the LOBSTER biogeochemical model. It includes single Phytoplankton and Zooplankton tracers which track the nitrogen (mmol N / m³).

The phytoplankton evolve like ∂ₜP = (1-γ)μP - Gₚ - mP² where μ is the growth rate (1/s) which is of the form phytoplankton_maximum_growth_rate × light limitation × nutrient limitation. The nutrient_limitation term can be extended to account for different nutrient formulations. The (1-γ) term is the assimilated fraction. Gₚ accounts for grazing by the zooplankton, and m is the quadratic mortality rate (1/s/mmol N/m³).

The zooplankton evolves like ∂ₜZ = αG - mZ² - μZ where G is the total grazing of phytoplankton and detritus with α the fraction assimilated, m is the quadratic mortality rate, and μ is the linear mortality rate taken to represent the excretion rate.

TwoParticleAndDissolved

TwoParticleAndDissolved defines the default living component for the LOBSTER biogeochemical model. It includes small and large particulate organic matter (sPOM and bPOM), and dissolved organic matter (DOM).

POM is produced by "waste" from the biological system (e.g. excretion from zooplankton or phytoplankton death) and is portioned into sPOM and bPOM as by the small_solid_waste_fraction parameter. sPOM is depleated by grazing by the biology (typically the zooplankton), and is remineralised. Calcite produced by the biology (from phytoplankton mortality) is assumed to all end up in bPOM.

DOM is produced by the biological system (e.g. excretion) and the breakdown of POM, and is broken down at a constant remineralisation rate.

Note: we refer to the particle classes as "small" and "large", but label them sPOM and bPOM for "small" and "big" as using lPOM may be confused for labile POM. We agree that this is a clumsy solution and would be happy to find an alternative.

TwoParticleAndDissolved(grid; 
                        small_particle_sinking_speed = 3.47e-5, # m/s
                        large_particle_sinking_speed = 200/day, # m/s
                        open_bottom = true,
                        kwargs...)

Construct an instance of TwoParticleAndDissolved specifying the small_particle_sinking_speed and large_particle_sinking_speed. If open_bottom is true then particles sink out of the bottom, if false then the sinking velocity smoothly goes to zero at the bottom to prevent the tracers leaving the domain.

VariableRedfieldDetritus

VariableRedfieldDetritus defines builds on the default TwoParticleAndDissolved component for the LOBSTER biogeochemical model, splitting each component into separate carbon and nitrogen compartements: sPOC, bPOC, DOC, sPON, bPON, and DON.

Standing alone this will produce identical behavior to the TwoParticleAndDissolved formulation, but if particles are injected with a different C:N ratio (for example) from the SugarKelp model, then the ratio will evolve differently.

VariableRedfieldDetritus(grid; 
                         small_particle_sinking_speed = 3.47e-5, # m/s
                         large_particle_sinking_speed = 200/day, # m/s
                         open_bottom = true,
                         kwargs...)

Construct an instance of TwoParticleAndDissolved specifying the small_particle_sinking_speed and large_particle_sinking_speed. If open_bottom is true then particles sink out of the bottom, if false then the sinking velocity smoothly goes to zero at the bottom to prevent the tracers leaving the domain.

Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) model.

This is not currently an official version supported by the PISCES community and is not yet verified to be capable of producing results mathcing that of the operational PISCES configuration. This is a work in progress, please open an issue or discusison if you'd like to know more.

Notes to developers

Part of the vision for this implementation of PISCES is to harness the features of Julia that would allow it to be fully modular. An obvious step to improve the ease of this would be to do some minor refactoring to group the phytoplankton classes, and zooplankton classes together, and for the other groups to generically call the whole lot. This may cause some issues with argument passing, and although it may not be the best way todo it my first thought is to pass them round as named tuples built from something like,

phytoplankton_tracers = phytoplankton_arguments(bgc.phytoplankton, args...)

DepthDependantSinkingSpeed(; minimum_speed = 30/day,
                             maximum_speed = 200/day,
                             maximum_depth = 500)

Returns sinking speed for particles which sink at minimum_speed in the surface ocean (the deepest of the mixed and euphotic layers), and accelerate to maximum_speed below that depth and maximum_depth.

ModelLatitude

Returns the latitude specified by the model grid (y).

PISCES(; grid::AbstractGrid{FT},
         phytoplankton = MixedMondoNanoAndDiatoms(),
         zooplankton = MicroAndMesoZooplankton(),
         dissolved_organic_matter = DissolvedOrganicCarbon(),
         particulate_organic_matter = TwoCompartementCarbonIronParticles(),
         
         nitrogen = NitrateAmmonia(),
         iron = SimpleIron(),
         silicate = Silicate(),
         oxygen = Oxygen(),
         phosphate = Phosphate(),
         
         inorganic_carbon = InorganicCarbon(),

         # from Aumount 2005 rather than 2015 since it doesn't work the other way around
         first_anoxia_thresehold = 6.0,
         second_anoxia_thresehold = 1.0,

         nitrogen_redfield_ratio = 16/122,
         phosphate_redfield_ratio = 1/122,
         
         mixed_layer_shear = 1.0,
         background_shear = 0.01, 
         
         latitude = PrescribedLatitude(45),
         day_length = day_length_function,
         
         mixed_layer_depth = Field{Center, Center, Nothing}(grid),
         euphotic_depth = Field{Center, Center, Nothing}(grid),

         silicate_climatology = ConstantField(7.5),

         mean_mixed_layer_vertical_diffusivity = Field{Center, Center, Nothing}(grid),
         mean_mixed_layer_light = Field{Center, Center, Nothing}(grid),

         carbon_chemistry = CarbonChemistry(),
         calcite_saturation = CenterField(grid),

         surface_photosynthetically_active_radiation = default_surface_PAR,

         light_attenuation =
           MultiBandPhotosyntheticallyActiveRadiation(; grid, 
                                                        surface_PAR = surface_photosynthetically_active_radiation),

         sinking_speeds = (POC = 2/day, 
                           # might be more efficient to just precompute this
                           GOC = Field(KernelFunctionOperation{Center, Center, Face}(DepthDependantSinkingSpeed(), 
                                                                                     grid, 
                                                                                     mixed_layer_depth, 
                                                                                     euphotic_depth))),
         open_bottom = true,

         scale_negatives = false,
         invalid_fill_value = NaN,
         
         sediment = nothing,
         particles = nothing,
         modifiers = nothing)

Constructs an instance of the PISCES biogeochemical model.

Keyword Arguments

• grid: (required) the geometry to build the model on
• phytoplankton: phytoplankton evolution parameterisation, defaults to nanophyto and diatom size classes with MixedMondo growth
• zooplankton: zooplankton evolution parameterisation, defaults to two class Z and M
• dissolved_organic_matter: parameterisaion for the evolution of dissolved organic matter (DOC)
• particulate_organic_matter: parameterisation for the evolution of particulate organic matter (POC, GOC, SFe, BFe, PSi, CaCO₃)
• nitrogen: parameterisation for the nitrogen compartements (NH₄ and NO₃)
• iron: parameterisation for iron (Fe), currently the "complex chemistry" of Aumount 2015 is not implemented
• silicate: parameterisaion for silicate (Si)
• oxygen: parameterisaion for oxygen (O₂)
• phosphate: parameterisaion for phosphate (PO₄)
• inorganic_carbon: parameterisation for the evolution of the inorganic carbon system (DIC and Alkalinity)
• first_anoxia_thresehold and second_anoxia_thresehold: thresholds in anoxia parameterisation
• nitrogen_redfield_ratio and phosphate_redfield_ratio: the assumed element ratios N/C and P/C
• mixed_layer_shear and background_shear: the mixed layer and background shear rates, TODO: move this to a computed field
• latitude: model latitude, should be PrescribedLatitude for RectilinearGrids and ModelLatitude for grids providing their own latitude
• day_length: parameterisation for day length based on time of year and latitude, you may wish to change this to (φ, t) -> 1day if you want to ignore the effect of day length, or something else if you're modelling a differen planet
• mixed_layer_depth: an AbstractField containing the mixed layer depth (to be computed during update state)
• euphotic: an AbstractField containing the euphotic depth, the depth where light reduces to 1/1000 of the surface value (computed during update state)
• silicate_climatology: an AbstractField containing the silicate climatology which effects the diatoms silicate half saturation constant
• mean_mixed_layer_vertical_diffusivity: an AbstractField containing the mean mixed layer vertical diffusivity (to be computed during update state)
• mean_mixed_layer_light: an AbstractField containing the mean mixed layer light (computed during update state)
• carbon_chemistry: the CarbonChemistry model used to compute the calicte saturation
• calcite_saturation: an AbstractField containing the calcite saturation (computed during update state)
• surface_photosynthetically_active_radiation: funciton for the photosynthetically available radiation at the surface
• light_attenuation: light attenuation model which integrated the attenuation of available light
• sinking_speed: named tuple of constant sinking speeds, or fields (i.e. ZFaceField(...)) for any tracers which sink (convention is that a sinking speed is positive, but a field will need to follow the usual down being negative)
• open_bottom: should the sinking velocity be smoothly brought to zero at the bottom to prevent the tracers leaving the domain
• scale_negatives: scale negative tracers?
• particles: slot for BiogeochemicalParticles
• modifiers: slot for components which modify the biogeochemistry when the tendencies have been calculated or when the state is updated

All parameterisations default to the operaitonal version of PISCES as close as possible.

Notes

Currently only MixedMondoPhytoplankton are implemented, and some work should be done to generalise the classes to a single phytoplankton if more classes are required (see OceanBioME.Models.PISCESModel docstring). Similarly, if a more generic particulate_organic_matter was desired a way to specify arbitary tracers for arguments would be required.

PrescribedLatitude

Returns the prescribed latitude rather than the model grid y position.

DissolvedOrganicCarbon

Parameterisation of dissolved organic matter which depends on a bacterial concentration.

TwoCompartementCarbonIronParticles

A quota parameterisation for particulate organic matter with two size classes, each with carbon and iron compartements, and a silicate compartement for the large size class.

Confusingly we decided to name these compartmenets POC and GOC for the small and large carbon classes, SFe and BFe for the small and ̶l̶a̶r̶g̶e̶ big iron compartements, and PSi for the ̶l̶a̶r̶g̶e̶ particulate silicon (not the phytoplankton silicon).

SimpleIron(; excess_scavenging_enhancement = 1000)

Parameterisation for iron evolution, not the "complex chemistry" model of Aumount et al, 2015. Iron is scavenged (i.e. perminemtly removed from the model) when the free iron concentration exeeds the ligand concentration at a rate modified by excess_scavenging_enhancement.

InorganicCarbon

Default parameterisation for DICandAlk`alinity evolution.

QualityDependantZooplankton

The PISCES zooplankton growth model where each class has preferences for grazing on nanophytoplankton (P), diatoms (D), microzooplankton (Z), and particulate organic matter (POC), and can flux feed on sinking particulates (POC and GOC).

This model assumes a fixed ratio for all other elements (i.e. N, P, Fe).

NutrientLimitedProduction

BaseProduction with light limitation moderated by nutrient availability. This is the "new production" PISCES phytoplankton growth rate model. Growth rate is of the form:

\[μ = μ₁f₁(τ)f₂(zₘₓₗ)(1-exp(-α θᶜʰˡ PAR / τ (bᵣ + μᵣ))) L.\]

Keyword Arguments

• base_growth_rate: the base growth rate, μ₀, in (1/s)
• temperatrue_sensetivity: temperature sensetivity parameter, b, giving μ₁ = μ₀ bᵀ where T is temperature
• dark_tollerance: the time that the phytoplankton survives in darkness below the euphotic layer, τᵈ (s)
• initial_slope_of_PI_curve: the relationship between photosynthesis and irradiance, α₀ (1/W/m²)
• low_light_adaptation: factor increasing the sensetivity of photosynthesis to irradiance, β, giving α = α₀(1 + exp(-PAR)), typically set to zero
• basal_respiration_rate: reference respiration rate, bᵣ (1/s)
• reference_growth_rate: reference growth rate, μᵣ (1/s)

MixedMondo

Holds the parameters for the PISCES mixed mondo phytoplankton parameterisation where nutrient limitation is modelled using the mondo approach for nitrate (NO₃), ammonia (NH₄), phosphate (PO₄), and silicate (Si), but the quota approach is used for iron (Fe) and light (PAR).

Therefore each class has a carbon compartement (generically I), chlorophyll (IChl), and iron (IFe), and may also have silicate (ISi) if the nutrient_limitation specifies that the growth is silicate limited, despite the fact that the silicate still limits the growth in a mondo fashion.

The growth_rate may be different parameterisations, currently either NutrientLimitedProduction or GrowthRespirationLimitedProduction, which represent the typical and newprod versions of PISCES.

NitrogenIronPhosphateSilicateLimitation

Holds the parameters for growth limitation by nitrogen (NO₃ and NH₄), iron (Fe), phosphate PO₄, and (optionally) silicate (Si) availability.

Silicate limitation may be turned off (e.g. for nanophytoplankton) by setting silicate_limited=false.

NutrientLimitedProduction

BaseProduction with light limitation moderated by nutrient availability. This is the "origional" PISCES phytoplankton growth rate model. Growth rate is of the form:

\[μ = μ₁f₁(τᵈ)f₂(zₘₓₗ)(1-exp(-α θᶜʰˡ PAR / τ μ₀ L)) L.\]

Keyword Arguments

• base_growth_rate: the base growth rate, μ₀, in (1/s)
• temperatrue_sensetivity: temperature sensetivity parameter, b, giving μ₁ = μ₀ bᵀ where T is temperature
• dark_tollerance: the time that the phytoplankton survives in darkness below the euphotic layer, τᵈ (s)
• initial_slope_of_PI_curve: the relationship between photosynthesis and irradiance, α₀ (1/W/m²)
• low_light_adaptation: factor increasing the sensetivity of photosynthesis to irradiance, β, giving α = α₀(1 + exp(-PAR)), typically set to zero

NitrateAmmonia

A parameterisation for the evolution of nitrate (NO₃) and ammonia (NH₄) where ammonia can be nitrified into nitrate, nitrate and ammonia are supplied by the bacterial degredation of dissolved organic matter, and consumed by phytoplankton. Additionally waste produces ammonia through various means.

Sugar kelp model of Broch and Slagstad (2012) and updated by Broch et al. (2013), Fossberg et al. (2018), and Broch et al. (2019).

Prognostic properties

• Area: A (dm²)
• Nitrogen reserve: N (gN/gSW)
• Carbon reserve: C (gC/gSW)

Tracer dependencies

• Nitrates: NO₃ (mmol N/m³)
• Ammonia: NH₄ (mmol N/m³)
• Photosynthetically available radiation: PAR (einstein/m²/day)
• Temperature: T (°C)

SugarKelp

Defines the parameters for SugarKelp biogeochemistry.

SugarKelpParticles(n; grid, kelp_parameters = NamedTuple(), kwargs...)

Sets up n sugar kelp BiogeochemicalParticles with default parameters except those specified in kelp_parameters. kwagrs are passed onto BiogeochemicalParticles.

CarbonChemistryModel to solve chemical equilibrium parameterisations

CarbonChemistry(FT = Float64; 
                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,
                first_virial_coefficient = PolynomialVirialCoefficientForCarbonDioxide(),
                cross_viral_coefficient = CrossVirialCoefficientForCarbonDioxide(),
                solver = NewtonRaphsonSolver())

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

julia> using OceanBioME

julia> carbon_chemistry = CarbonChemistry()
`CarbonChemistry` model which solves for pCO₂ and pH

julia> pCO₂ = carbon_chemistry(; DIC = 2000.0, Alk = 2000.0, T = 10.0, S = 35.0)
1308.1474527899106

julia> pH = carbon_chemistry(; DIC = 2000.0, Alk = 2000.0, T = 10.0, S = 35.0, output = Val(:pHᶠ))
7.502532746463654

julia> pCO₂_higher_pH = carbon_chemistry(; DIC = 2000.0, T = 10.0, S = 35.0, pH = 7.5)
1315.7136384737507

(p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing,
                       output = Val(:fCO₂),
                       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,
                       initial_pH_guess = 8)

Calculates fCO₂ in sea water with DIC, Alkalinity, Temperature, and Salinity unless pH is specified, in which case intermediate computation of pH is skipped and pCO₂ is calculated from the DIC, T, S and pH.

DIC is expected in mmol C/m³, Alk meq/m³, T in °C, and S in PSU.

When pH is specified the free pH (i.e. -log[H⁺]) is expected.

Alternativly pCO₂, and free, total, or sea water pH may be returned by setting output to Val(:pCO₂), Val(:pHᶠ), Val(:pHᵗ), or Val(:pHˢ), which will return X in Val(:X) instead of fCO₂.

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

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

K1(FT = Float64;
   constant =  61.2172,
   inverse_T = -3633.86,
   log_T = -9.67770,
   S =  0.011555,
   S² = -0.0001152,
   pressure_correction = PressureCorrection(FT; 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).

K2(FT = Float64;
   constant = -25.9290,
   inverse_T = -471.78,
   log_T = 3.16967,
   S = 0.01781,
   S² = -0.0001122,
   pressure_correction = PressureCorrection(FT; 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).

KB(FT = Float64;
   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(FT; 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).

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

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.

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

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

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

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

PressureCorrection(FT=Float64;
                   a₀, a₁, a₂,
                   b₀, b₁, b₂,
                   R = 83.14472)

Parameterisation for the pressure effect on thermodynamic constants.

Form from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341).

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

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

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

KSParagonite(; thermconstant = -171.945, thermT = -0.077993, therminverseT = 2903.293, thermlogT = 71.595, seasqrtS = -0.068393, seaTsqrtS = 0.0017276, seainverseTsqrtS = 88.135, seaS = -0.10018, seaSsqrtS³ = 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_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).

Light attenuation by chlorophyll as described by Karleskind et al. (2011) (implemented as TwoBand) and Morel (1988) (as MultiBand).

MultiBandPhotosyntheticallyActiveRadiation{F, FN, K, E, C, SPAR, SPARD}

Light attenuation model with multiple wave length bands where each band (i) is attenuated like:

∂PARᵢ(z)/∂z = PARᵢ(kʷ(i) + χ(i)Chl(z)ᵉ⁽ⁱ⁾)

Where kʷ(i) is the band specific water attenuation coefficient, e(i) the chlorophyll exponent, and χ(i) the chlorophyll attenuation coefficient.

When the fields are called with biogeochemical_auxiliary_fields an additional field named PAR is also returned which is a sum of the bands.

MultiBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, 
                                             bands = ((400, 500), (500, 600), (600, 700)), #nm
                                             base_bands = MOREL_λ,
                                             base_water_attenuation_coefficient = MOREL_kʷ,
                                             base_chlorophyll_exponent = MOREL_e,
                                             base_chlorophyll_attenuation_coefficient = MOREL_χ,
                                             field_names = [par_symbol(n, length(bands)) for n in 1:length(bands)],
                                             surface_PAR = default_surface_PAR)

Returns a MultiBandPhotosyntheticallyActiveRadiation attenuation model of surface_PAR in divided into bands by surface_PAR_division.

The attenuation morelcoefficients are computed from `basewaterattenuationcoefficient,basechlorophyllexponent, andbasechlorophyllattenuationcoefficientwhich should be arrays of the coefficients atbasebands` wavelengths.

The returned field_names default to PAR₁, PAR₂, etc., but may be specified by the user instead.

Keyword Arguments

• grid: grid for building the model on
• water_red_attenuation, ..., phytoplankton_chlorophyll_ratio: parameter values
• surface_PAR: function (or array in the future) for the photosynthetically available radiation at the surface, which should be f(x, y, t) where x and y are the native coordinates (i.e. meters for rectilinear grids and latitude/longitude as appropriate)

PrescribedPhotosyntheticallyActiveRadiation(fields)

PrescribedPhotosyntheticallyActiveRadiation returns "prescribed" PAR fields which are user specified, e.g. they may be FunctionFields or ConstantFields.

fields may either be an AbstractField or a NamedTuple of names and fields which will be returned in biogeochemical_auxiliary_fields, if only one field is present the field will be named PAR.

TwoBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, 
                                           water_red_attenuation::FT = 0.225, # 1/m
                                           water_blue_attenuation::FT = 0.0232, # 1/m
                                           chlorophyll_red_attenuation::FT = 0.037, # 1/(m * (mgChl/m³) ^ eʳ)
                                           chlorophyll_blue_attenuation::FT = 0.074, # 1/(m * (mgChl/m³) ^ eᵇ)
                                           chlorophyll_red_exponent::FT = 0.629,
                                           chlorophyll_blue_exponent::FT = 0.674,
                                           pigment_ratio::FT = 0.7,
                                           phytoplankton_chlorophyll_ratio::FT = 1.31,
                                           surface_PAR::SPAR = default_surface_PAR)

Keyword Arguments

• grid: grid for building the model on
• water_red_attenuation, ..., phytoplankton_chlorophyll_ratio: parameter values
• surface_PAR: function (or array in the future) for the photosynthetically available radiation at the surface, which should be f(x, y, t) where x and y are the native coordinates (i.e. meters for rectilinear grids and latitude/longitude as appropriate)

InstantRemineralisation

Hold the parameters and fields a simple sediment model where sinking organic carbon is "instantly remineralised" and either returned to the domain as remineralisation_reciever (typically NH₄), or permanently stored in the sediment.

The "burial efficiency" (the fraction permanently stored) is from Dunne et al. (2007), and varies with the sinking flux.

struct SimpleMultiG

Hold the parameters and fields for a simple "multi G" single-layer sediment model. Based on the Level 3 model described by Soetaert et al. (2000).

InstantRemineralisationSediment(grid;
                                sinking_tracers = (:P, :D), 
                                remineralisation_reciever = :N,
                                burial_efficiency_constant1 = 0.013,
                                burial_efficiency_constant2 = 0.53,
                                burial_efficiency_half_saturation = 7.0 / 6.56,
                                kwargs...)

Return a single-layer instant remineralisation sediment model where the sinking_tracers are instantly remineralised and returned to remineralisation_reciever with a small fraction permanently buried with efficiency:

e = a + b * f / (k + f)²

where a is burial_efficiency_constant1, b is burial_efficiency_constant2, and k is the burial_efficiency_half_saturation.

kwargs... are BiogeochemicalSediment key word arguments.

Example

using OceanBioME, Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))

sediment_model = InstantRemineralisationSediment(grid)

biogeochemistry = NPZD(; grid, sediment_model)

using OceanBioME, Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))

sediment_model = InstantRemineralisationSediment(grid; 
                                                 sinking_tracers = (:sPOM, :bPOM),
                                                 remineralisation_reciever = :NH₄)

biogeochemistry = LOBSTER(; grid, sediment_model)

SimpleMultiGSediment(grid;
                     fast_decay_rate = 2/day,
                     slow_decay_rate = 0.2/day,
                     fast_redfield = 0.1509,
                     slow_redfield = 0.13,
                     fast_fraction = 0.74,
                     slow_fraction = 0.26,
                     refactory_fraction = 0.1,
                     sedimentation_rate = 982 * abs(znode(1, 1, 1, grid, Center(), Center(), Center())) ^ (-1.548), # cm/year, incorrect for D < 100m
                     anoxia_half_saturation = 1.0, # mmol/m³ (arbitarily low)
                     nitrate_oxidation_params = on_architecture(architecture(grid), (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232)),
                     denitrification_params = on_architecture(architecture(grid), (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662)),
                     anoxic_params = on_architecture(architecture(grid), (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143)),
                     solid_dep_params = on_architecture(architecture(grid), (0.233, 0.336, 982.0, - 1.548)),
                     sinking_nitrogen = (:sPOM, :bPOM),
                     sinking_carbon = nothing,
                     sinking_redfield = ifelse(isnothing(sinking_carbon), convert(eltype(grid), 6.56), nothing),
                     kwargs...)

Return a single-layer "multi G" sediment model (SimpleMultiG) on grid, where parameters can be optionally specified.

The model is a single layer (i.e. does not include porous diffusion) model with three classes of sediment organic matter which decay at three different rates (fast, slow, refactory). The nitrification/denitrification/anoxic mineralisation fractions default to the parameterisation of Soetaert et al. 2000; doi:10.1016/S0012-8252(00)00004-0.

This model has not yet been validated or compared to observational data. The variety of degridation processes is likely to be strongly dependent on oxygen availability (see https://bg.copernicus.org/articles/6/1273/2009/bg-6-1273-2009.pdf) so it will therefore be important to also thoroughly validate the oxygen model (also currently limited).

Example

julia> using OceanBioME, Oceananigans

julia> grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));

julia> sediment = SimpleMultiGSediment(grid)
`BiogeochemicalSediment` with `Single-layer multi-G sediment model (Float64)` biogeochemsitry
    Prognostic fields: (:Ns, :Nf, :Nr)
    Tracked fields: (:NO₃, :NH₄, :O₂, :sPOM, :bPOM)
    Coupled fields: (:NO₃, :NH₄, :O₂)

julia> biogeochemistry = LOBSTER(; grid, sediment, detritus = TwoParticleAndDissolved(grid; open_bottom=true))
LOBSTER model (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: `BiogeochemicalSediment` with `Single-layer multi-G sediment model (Float64)` biogeochemsitry
 Particles: Nothing
 Modifiers: Nothing

GasExchangeModel to solve chemical equilibrium parameterisations

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.

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.

PartiallySolubleGas(; air_concentration, solubility)

Parameterises the available concentration of a gas dissolving in water in the form $\alpha C_a$ where $lpha$ is the Ostwald solubility coeffieient and $C_a$ is the concentration in the air.

Wanninkhof92Solubility

Parameterises the Ostwald solubility coefficient as given in Wanninkhof, 1992.

CarbonDioxideGasExchangeBoundaryCondition(FT = Float64; 
                                          carbon_chemistry = CarbonChemistry(FT),
                                          transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = CarbonDioxidePolynomialSchmidtNumber(FT)),
                                          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 Tupleof symbols specifying the name of the silicate and phosphate tracers, or aNamedTupleof values for thecarbon_chemistry` model.

kwargs are passed on to GasExchangeBoundaryCondition.

Note: The model always requires T, S, DIC, and Alk to be present in the model.

CarbonDioxidePolynomialSchmidtNumber(FT = Float64; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555)

Schmidt number parameterisation Wanninkhof, 2014 for sea water

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

OxygenGasExchangeBoundaryCondition(FT = Float64; 
                                   transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = OxygenPolynomialSchmidtNumber(FT)),
                                   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.

OxygenPolynomialSchmidtNumber(FT = Float64; a = 1953.4, b = - 128.0, c = 3.9918, d = -0.050091)

Schmidt number parameterisation Wanninkhof, 2014 for sea water

SchmidtScaledTransferVelocity(; schmidt_number, 
                                base_transfer_velocity = Ho06())

Returns a model for gas transfer velocity which depends on the u₁₀, the 10m-wind, and Temperature. 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.

CCMP2(FT = Float64; 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).

ERA5(FT = Float64; 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).

Ho06(FT = Float64; 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.

JRA55(FT = Float64; 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).

McGillis01(FT = Float64; 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.

NCEP1(FT = Float64; 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).

Nightingale00(FT = Float64; 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 (?).

Sweeny07(FT = Float64; scale_factor = 0.27 / hour / 100)

Quadratic k₆₆₀ parameterisation of Sweeny et al. (2007) suitable for the NCEP/NCAR reanalysis 1 product

Wanninkhof09(FT = Float64; 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

Wanninkhof14(FT = Float64; scale_factor = 0.251 / hour / 100)

Quadratic k₆₆₀ parameterisation of Wanninkhof et al (2014) suitable for the Cross-Calibrated Multi-Platform (CCMP) Winds product

Wanninkhof99(FT = Float64; scale_factor = 0.0283 / hour / 100)

Cubic k₆₆₀ parameterisation of Wanninkhof & McGillis (1999) suitable for short term, in situ wind products.

Integrate biogeochemical models on a single point

BoxModel(; biogeochemistry,
           forcing = NamedTuple(),
           timestepper = :RungeKutta3,
           clock = Clock(; time = 0.0),
           prescribed_tracers::PT = NamedTuple())

Constructs a box model of a biogeochemistry model. Once this has been constructed you can set initial condiitons by set!(model, X=1.0...).

Keyword Arguments

• biogeochemistry: (required) an OceanBioME biogeochemical model, most models must be passed a grid which can be set to a BoxModelGrid for box models
• 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: named tuple of tracer names and function (f(t)) prescribing tracer values

set!(model::BoxModel; kwargs...)

Set the values for a BoxModel

Arguments

• model - the model to set the arguments for

Keyword Arguments

• variables and value pairs to set

Store previous source terms before updating them.

BiogeochemicalParticles(number;
                        grid,
                        biogeochemistry,
                        advection = LagrangianAdvection(),
                        timestepper = ForwardEuler,
                        field_interpolation = NearestPoint(),
                        scalefactors = ones(number))

Creates number particles with biogeochemistry on grid, advected by advection which defaults to LagrangianAdvection (i.e. they comove with the water). The biogeochemistry is stepped by timestepper and tracer fields are interpolated by field_interpolation, which defaults to directly reading the nearest center point and taking up/depositing in the same.

Particles can also have a scalefactor which scales their tracer interaction (e.g. to mimic the particle representing multiple particles).

ForwardEuler

Step particle biogeochemistry with a ForwardEuler methods with Δt from the physical model substep.

LagrangianAdvection

Specifies that particles should move in a purley lagrangian mannor.

NearestPoint

Specifies that tracer values should be taken from the nearst center point.