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.

NutrientsPlanktonDetritus(grid;
                          nutrients = Nutrients(nothing, nothing, nothing, nothing),
                          plankton = Abiotic(),
                          detritus = InstantRemineralisationDetritus(),
                          inorganic_carbon = nothing,
                          oxygen = nothing,
                          light_attenuation = nothing,
                          sediment = nothing,
                          scale_negatives = false,
                          invalid_fill_value = NaN,
                          particles = nothing,
                          modifiers = nothing)

Construct a biogeochemical model in the modular Nutrients-Plankton-Detritus (NPD) framework by assembling it from pluggable components. Each slot may be swapped independently, and the preset LOBSTER, NPZD, and ImplicitBiology constructors are just particular choices of these components.

The set of tracers the model evolves is determined by the components you choose.

Keyword Arguments

• grid: (required) the geometry to build the model on, required to configure sinking speeds
• nutrients: the inorganic nutrient pool(s) that limit growth; a Nutrients grouping of nitrogen, phosphate, iron, and silicate components. Defaults to no explicit nutrients (all slots nothing), in which case nutrients are implicitly conserved and not tracked
• plankton: the planktonic (living) component, defaults to Abiotic (no biology). Options include PhytoZoo and ImplicitProductivity
• detritus: the non-living organic component, defaults to InstantRemineralisationDetritus which returns waste straight to the nutrient pool. Options include Detritus, DissolvedParticulate, and CarbonNitrogenDissolvedParticulate
• inorganic_carbon: optional inorganic carbon component, defaults to nothing; can be a CarbonateSystem which adds dissolved inorganic carbon (DIC) and alkalinity (Alk)
• oxygen: optional oxygen component, defaults to nothing; can be an Oxygen which adds oxygen (O₂)
• light_attenuation: light attenuation model which integrates the attenuation of available light
• sediment: slot for a sediment model (AbstractSediment)
• scale_negatives: whether to add a ScaleNegativeTracers modifier to keep tracers non-negative
• invalid_fill_value: the value used to fill invalid tracer values when scale_negatives is true
• particles: slot for BiogeochemicalParticles
• modifiers: slot for components which modify the biogeochemistry after the tendencies have been calculated or when the state is updated

ImplicitBiology(grid; limiting_nutrients = (:nitrate, :iron, :phosphate), open_bottom = true, kwargs...)

Construct the ImplicitBiology preset of the NutrientsPlanktonDetritus framework: a model which computes community productivity limited by available nutrients and light without explicitly tracking planktonic biomass (see ImplicitProductivity). It is suited to large-scale or long-timescale simulations where resolving full plankton dynamics is too expensive.

By default it couples the Nutrients selected by limiting_nutrients, an ImplicitProductivity plankton, a two-class DissolvedParticulate detritus (DOP/POP), a CarbonateSystem, and a PrescribedAttenuationPAR light model.

Keyword Arguments

• grid: (required) the geometry to build the model on
• limiting_nutrients: a tuple of the nutrients that limit productivity, drawn from :nitrate (or :ammonia to split nitrogen into NO₃/NH₄), :phosphate, and :iron
• open_bottom: whether detritus can sink out of the bottom of the domain
• surface_PAR: the surface photosynthetically active radiation passed to the default light model
• any other keyword argument is forwarded to NutrientsPlanktonDetritus (e.g. oxygen, nutrients, plankton, detritus, light_attenuation)

LOBSTER(grid; limiting_nutrients = (:nitrate, :ammonia), open_bottom = true, kwargs...)

Construct the LOBSTER preset of the NutrientsPlanktonDetritus framework. By default it splits nitrogen into nitrate (NO₃) and ammonia (NH₄) via NitrateAmmonia, couples a PhytoZoo plankton (phytoplankton P and zooplankton Z), a dissolved-and-two- particulate DissolvedParticulate detritus (DOM, sPOM, bPOM), and the default two-band light model. Pass inorganic_carbon = CarbonateSystem() and/or oxygen = Oxygen() to add carbonate chemistry and oxygen.

Keyword Arguments

• grid: (required) the geometry to build the model on
• limiting_nutrients: a tuple of the nutrients that limit growth, drawn from :nitrate, :ammonia, :phosphate, and :iron
• open_bottom: whether detritus can sink out of the bottom of the domain
• surface_PAR: the surface photosynthetically active radiation passed to the default light model
• any other keyword argument is forwarded to NutrientsPlanktonDetritus (e.g. inorganic_carbon, oxygen)

NPZD(grid; limiting_nutrients = (:nitrate,), open_bottom = true, kwargs...)

Construct the NPZD (Nutrient-Phytoplankton-Zooplankton-Detritus) preset of the NutrientsPlanktonDetritus framework. It couples the Nutrients selected by limiting_nutrients with a PhytoZoo plankton (phytoplankton P and zooplankton Z) parameterised after Kuhn et al. (2015), a single-class Detritus pool (D), and the default two-band light model.

Keyword Arguments

• grid: (required) the geometry to build the model on
• limiting_nutrients: a tuple of the nutrients that limit growth, drawn from :nitrate (or :ammonia to split nitrogen into NO₃/NH₄), :phosphate, and :iron
• open_bottom: whether detritus can sink out of the bottom of the domain
• surface_PAR: the surface photosynthetically active radiation passed to the default light model
• any other keyword argument is forwarded to NutrientsPlanktonDetritus

NitrateAmmonia(; nitrification_rate = 5.8e-7)

A nitrogen component for the nitrogen slot of Nutrients that resolves nitrate (NO₃) and ammonia (NH₄) as separate tracers, with ammonia converted to nitrate by nitrification. Plankton may then take up nitrate and ammonia independently and preferentially.

Keyword Arguments

• nitrification_rate: the rate at which ammonia is nitrified to nitrate (1/s)

Nutrients(; nitrogen = nothing, phosphate = nothing, iron = nothing, silicate = nothing)

Group the inorganic nutrient pools that limit plankton growth in a NutrientsPlanktonDetritus model. Each of the four slots may be nothing (the nutrient is not tracked and is implicitly conserved), a SingleTracerNutrient (N, PO₄, Fe, or Si), or a more detailed component such as NitrateAmmonia for the nitrogen slot.

The tracers added to the model are the union of those required by the non-nothing slots.

Keyword Arguments

• nitrogen: the nitrogen pool, e.g. N (single tracer) or NitrateAmmonia (NO₃/NH₄)
• phosphate: the phosphate pool, e.g. PO₄
• iron: the iron pool, e.g. Fe
• silicate: the silicate pool, e.g. Si

SingleTracerNutrient

An @enum of the single-tracer inorganic nutrients that can fill a Nutrients slot: N (nitrogen), PO₄ (phosphate), Fe (iron), and Si (silicate). Each introduces one tracer of the same name and is the simplest way to represent a limiting nutrient.

Abiotic()

An empty plankton component for the plankton slot of a NutrientsPlanktonDetritus model. It adds no tracers and produces no biological source or sink terms — all uptake and waste fluxes are zero. It is the default plankton and is useful for purely abiotic experiments (e.g. carbonate chemistry or gas exchange with no biology).

ImplicitProductivity

An implicit plankton which computes a community productivity limited by available nutrients and light, but not tracking any planktonic biomass.

Mainly from MITGCM https://web.mit.edu/globalchange/www/MITJPSPGC_Rpt122.pdf

Units: mmol P / m³ Tracers: none Default limiting nutrients: N or NO₃, P, Fe

PhytoZoo([FT = Float64;] kwargs...)
PhytoZoo(grid; phytoplankton_sinking_speed = 0, zooplankton_sinking_speed = 0, open_bottom = true, kwargs...)

A phytoplankton-and-zooplankton plankton component for the plankton slot of a NutrientsPlanktonDetritus model, as used by the LOBSTER and NPZD presets. It adds phytoplankton (P) and zooplankton (Z) tracers (plus a temperature tracer T when a temperature_coefficient is supplied). Phytoplankton grow with light- and nutrient-limited uptake, are grazed by zooplankton, and both pools produce dissolved/solid waste.

Use the grid method when the plankton should sink; it configures the sinking-velocity fields from phytoplankton_sinking_speed/zooplankton_sinking_speed and forwards the rest to the base constructor.

Keyword Arguments

• nutrient_half_saturations: a NamedTuple of half-saturation constants keyed by limiting nutrient (e.g. nitrate, ammonia, iron); its keys set which nutrients limit growth
• light_limitation: the light-limitation functional form (MondoLightLimitation() or AnalyticalLightLimitation()), with light_half_saturation (W/m²)
• phytoplankton_maximum_growth_rate, phytoplankton_exudation_fraction, phytoplankton_mortality_rate, phytoplankton_solid_waste_fraction: phytoplankton growth and loss parameters
• maximum_grazing_rate, grazing_half_saturation, preference_for_phytoplankton, zooplankton_assimilation_fraction, zooplankton_mortality_rate, zooplankton_excretion_rate: zooplankton grazing and loss parameters
• temperature_coefficient: optional Q₁₀ temperature dependence (off by default; adds a T tracer)
• carbon_ratio, iron_ratio, phosphate_ratio, rain_ratio, chlorophyll_ratio: elemental and chlorophyll ratios relative to nitrogen
• see the constructor definition for the full list of tunable parameters and their default values/units

CarbonNitrogenDissolvedParticulate(grid; kwargs...)

A variable-Redfield detritus component for the detritus slot of a NutrientsPlanktonDetritus model that tracks carbon and nitrogen separately in dissolved, small-particulate, and large-particulate classes. It adds the tracers DON, DOC, sPON, sPOC, bPON, and bPOC; the two particulate classes sink, and each pool remineralises to inorganic nutrients and carbon (partly via the dissolved pool).

Keyword Arguments

• grid: (required) the geometry, needed to configure the sinking-speed fields
• dissolved_remineralisation_rate, small_particle_remineralisation_rate, large_particle_remineralisation_rate: per-class remineralisation rates (1/s)
• small_fraction_of_solid_waste: the fraction of solid plankton waste routed to the small particulate class (the rest goes to the large class)
• small_particle_remineralisation_dissolved_fraction, large_particle_remineralisation_dissolved_fraction: the fraction of each particulate class's remineralisation that passes through the dissolved pool
• sinking_speeds: a NamedTuple (sPO = …, bPO = …) of the small/large particle sinking speeds (m/s)
• open_bottom: whether particulate detritus can sink out of the bottom of the domain

Detritus(grid; sinking_speed = 2.7489/day, open_bottom = true, remineralisation_rate = 0.1213/day)

A single-class detritus component (after Kuhn et al., 2015) for the detritus slot of a NutrientsPlanktonDetritus model. It adds one sinking detritus tracer D which accumulates plankton waste and grazing residue and is remineralised back to the nutrient pool.

Keyword Arguments

• grid: (required) the geometry, needed to configure the sinking-speed field
• sinking_speed: the downward sinking speed of detritus (m/s)
• open_bottom: whether detritus can sink out of the bottom of the domain
• remineralisation_rate: the rate at which detritus is remineralised to inorganic nutrients (1/s)

DissolvedParticulate(grid, dissolved_names = :DOM, particulate_names = (:sPOM, :bPOM); kwargs...)

A single-element detritus component for the detritus slot of a NutrientsPlanktonDetritus model that splits organic matter into any number of dissolved classes (which do not sink) and particulate classes (which sink). The tracers added are the dissolved_names and particulate_names (e.g. the default DOM, sPOM, bPOM). Plankton dissolved/solid waste is partitioned between the classes, particulate matter sinks and remineralises, and a fraction of particulate remineralisation may be routed through the dissolved pool.

Keyword Arguments

• grid: (required) the geometry, needed to configure the sinking-speed fields
• dissolved_names, particulate_names: the tracer name(s) of each dissolved / particulate class
• dissolved_remineralisation_rate, particulate_remineralisation_rate: per-class remineralisation rates (1/s)
• dissolved_waste_partitioning, particulate_waste_partitioning: how dissolved / solid plankton waste is split across the classes (each sums to 1)
• dissolved_fraction_of_remineralisation: the fraction of each particulate class's remineralisation that passes through the dissolved pool rather than directly to inorganic nutrients
• sinking_speeds: the downward sinking speed of each particulate class (m/s)
• open_bottom: whether particulate detritus can sink out of the bottom of the domain

InstantRemineralisationDetritus()

The simplest detritus component for the detritus slot of a NutrientsPlanktonDetritus model. It adds no tracers of its own: plankton dissolved and solid waste is instantaneously remineralised straight back to the inorganic nutrient pool (and, when present, inorganic carbon and calcite). This is the default detritus and is useful for closed-budget box models or when organic matter export is not of interest.

CarbonateSystem(replicates = 1)

The inorganic carbon component for the inorganic_carbon slot of a NutrientsPlanktonDetritus model. It adds dissolved inorganic carbon (DIC) and alkalinity (Alk) tracers whose evolution is driven by primary production, remineralisation of organic waste, and calcite production/dissolution (implicit calcite: calcite is not tracked as its own tracer).

Passing replicates > 1 manifests replicates independent copies of the carbonate tracers (DIC1, Alk1, DIC2, …), which is useful for ensemble or perturbation experiments; each replicate evolves with the same tendency as the base DIC/Alk.

Oxygen([FT = Float64;]
       production_oxygen_carbon_ratio = 131/122,
       nitrification_oxygen_carbon_ratio = 31/122)

An optional oxygen component for the oxygen slot of any NutrientsPlanktonDetritus model (including the LOBSTER, NPZD, and ImplicitBiology presets). It adds one tracer, oxygen (O₂), which is produced by photosynthesis and consumed by remineralisation of organic waste and by nitrification.

Oxygen concentration is only one-way coupled to the rest of the biogeochemistry and does not affect any other groups (e.g. low oxygen does not reduce zooplankton growth). To capture such an effect a different plankton component would be needed.

Keyword Arguments

• production_oxygen_carbon_ratio: moles of O₂ produced per mole of carbon fixed (mol O₂ / mol C)
• nitrification_oxygen_carbon_ratio: moles of O₂ consumed per mole of carbon during nitrification (mol O₂ / mol C)

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 matching that of the operational PISCES configuration. This is a work in progress, please open an issue or discussion 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 = TwoCompartmentCarbonIronParticles(),
         
         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_threshold = 6.0,
         second_anoxia_threshold = 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: parameterisation 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 compartments (NH₄ and NO₃)
• iron: parameterisation for iron (Fe), currently the "complex chemistry" of Aumount 2015 is not implemented
• silicate: parameterisation for silicate (Si)
• oxygen: parameterisation for oxygen (O₂)
• phosphate: parameterisation for phosphate (PO₄)
• inorganic_carbon: parameterisation for the evolution of the inorganic carbon system (DIC and Alk)
• first_anoxia_threshold and second_anoxia_threshold: 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 different 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 calcite saturation
• calcite_saturation: an AbstractField containing the calcite saturation (computed during update state)
• surface_photosynthetically_active_radiation: function for the photosynthetically available radiation at the surface
• light_attenuation: light attenuation model which integrated the attenuation of available light
• sinking_speeds: 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 operational 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 arbitrary 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.

TwoCompartmentCarbonIronParticles

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

Confusingly we decided to name these compartments POC and GOC for the small and large carbon classes, SFe and BFe for the small and ̶l̶a̶r̶g̶e̶ big iron compartments, 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. permanently removed from the model) when the free iron concentration exceeds the ligand concentration at a rate modified by excess_scavenging_enhancement.

InorganicCarbon

Default parameterisation for DICandAlk` 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)
• temperature_sensitivity: temperature sensitivity parameter, b, giving μ₁ = μ₀ bᵀ where T is temperature
• dark_tolerance: 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 sensitivity 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 compartment (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 "original" 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)
• temperature_sensitivity: temperature sensitivity parameter, b, giving μ₁ = μ₀ bᵀ where T is temperature
• dark_tolerance: 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 sensitivity 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 degradation 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.

Alternatively 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: the photosynthetically available radiation at the surface, by default, assumed to have the 'continuous form' condition(x, y t)

If parameters is not nothing, then surface_PAR has the form func(x, y, t, parameters).

If discrete_form = true, surface_PAR is assumed to have the "discrete form",

condition(i, j, grid, clock, fields)

where i, and j are indices that vary along the boundary. If discrete_form = true and parameters is not nothing, the function condition is called with

condition(i, j, grid, clock, fields, parameters)

PrescribedAttenuationPAR(grid, surface_PAR;
                         attenuation = 0.1,
                         attenuation_parameters = nothing,
                         attenuation_discrete_form = false,
                         surface_parameters = nothing,
                         surface_discrete_form = false)

Construct a light attenuation model (aliased PrescribedAttenuationPAR) which integrates a prescribed attenuation coefficient downward from the surface to give the photosynthetically active radiation (PAR), rather than computing attenuation from the chlorophyll concentration. It is the default light model for ImplicitBiology and is useful when no chlorophyll-based feedback on light is desired.

At depth $z$ the PAR is $\mathrm{PAR}(z) = \mathrm{PAR}_0 \exp\left(-\int_z^0 k \, dz'\right)$ where $k$ is the attenuation coefficient.

Arguments

• grid: the geometry to build the PAR field on
• surface_PAR: the photosynthetically active radiation at the surface; a number, or a function of the form f(x, y, t) (or the "discrete form" when surface_discrete_form = true)

Keyword Arguments

• attenuation: the attenuation coefficient $k$ (1/m); a number, or a function evaluated per cell
• attenuation_parameters, attenuation_discrete_form: parameters and form for attenuation when it is a function
• surface_parameters, surface_discrete_form: parameters and form for surface_PAR when it is a function

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 = 0.225, # 1/m
                                           water_blue_attenuation = 0.0232, # 1/m
                                           chlorophyll_red_attenuation = 0.037, # 1/(m * (mgChl/m³) ^ eʳ)
                                           chlorophyll_blue_attenuation = 0.074, # 1/(m * (mgChl/m³) ^ eᵇ)
                                           chlorophyll_red_exponent = 0.629,
                                           chlorophyll_blue_exponent = 0.674,
                                           pigment_ratio = 0.7,
                                           phytoplankton_chlorophyll_ratio = 1.31,
                                           surface_PAR = default_surface_PAR)

Keyword Arguments

• grid: grid for building the model on
• water_red_attenuation, ..., phytoplankton_chlorophyll_ratio: parameter values
• surface_PAR: the photosynthetically available radiation at the surface, by default, assumed to have the 'continuous form' condition(x, y t)

If parameters is not nothing, then surface_PAR has the form func(x, y, t, parameters).

If discrete_form = true, surface_PAR is assumed to have the "discrete form",

condition(i, j, grid, clock, fields)

where i, and j are indices that vary along the boundary. If discrete_form = true and parameters is not nothing, the function condition is called with

condition(i, j, grid, clock, fields, parameters)

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))
ERROR: UndefVarError: `TwoParticleAndDissolved` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] top-level scope
   @ none:1

GasExchangeModel to solve chemical equilibrium parameterisations

CarbonDioxideConcentration(FT = Float64;
                           carbon_chemistry::CC,
                           first_virial_coefficient::FV = PolynomialVirialCoefficientForCarbonDioxide{FT}(),
                           cross_virial_coefficient::CV = CrossVirialCoefficientForCarbonDioxide{FT}(),
                           air_pressure::AP = one(FT), # atm
                           silicate_and_phosphate_names::SP = nothing,
                           DIC = :DIC,
                           Alk = :Alk)

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.

DIC and Alk specify the tracer names of the DIC and alkalinity in the model.

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.