Library

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

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_model: 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

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(model::UnderlyingBiogeochemicalModel; warn = false)

Construct a modifier to scale the conserved tracers in model.

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.

conserved_tracers(model::UnderlyingBiogeochemicalModel)

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,
                                           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_model::LA =
                                               TwoBandPhotosyntheticallyActiveRadiation(; grid,
                                                                                          surface_PAR = surface_photosynthetically_active_radiation),
                                           sediment_model::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_model: light attenuation model which integrated the attenuation of available light
• sediment_model: slot for AbstractSediment
• 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.

Tracers

• Nitrates: NO₃ (mmol N/m³)
• Ammonia: NH₄ (mmol N/m³)
• Phytoplankton: P (mmol N/m³)
• Zooplankton: Z (mmol N/m³)
• Small (slow sinking) particulate organic matter: sPOM (mmol N/m³)
• Large (fast sinking) particulate organic matter: bPOM (mmol N/m³)
• Dissolved organic matter: DOM (mmol N/m³)

Optional tracers

Carbonate chemistry

• Dissolved inorganic carbon: DIC (mmol C/m³)
• Alkalinity: Alk (meq/m³)

Oxygen chemistry

• Oxygen: O₂ (mmol O₂/m³)

Variable redfield

• Small (slow sinking) particulate organic matter carbon content: sPOC (mmol C/m³)
• Large (fast sinking) particulate organic matter carbon content: bPOC (mmol C/m³)
• Dissolved organic matter carbon content: DOC (mmol C/m³)
• When this option is enabled then the usual sPOM and bPOM change to sPON and bPON as they explicitly represent the nitrogen contained in the particulate matter

Required submodels

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

For optional tracers:

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

LOBSTER(; grid,
          phytoplankton_preference::FT = 0.5,
          maximum_grazing_rate::FT = 9.26e-6, # 1/s
          grazing_half_saturation::FT = 1.0, # mmol N/m³
          light_half_saturation::FT = 33.0, # W/m² (?)
          nitrate_ammonia_inhibition::FT = 3.0,
          nitrate_half_saturation::FT = 0.7, # mmol N/m³
          ammonia_half_saturation::FT = 0.001, # mmol N/m³
          maximum_phytoplankton_growthrate::FT = 1.21e-5, # 1/s
          zooplankton_assimilation_fraction::FT = 0.7,
          zooplankton_mortality::FT = 2.31e-6, # 1/s/mmol N/m³
          zooplankton_excretion_rate::FT = 5.8e-7, # 1/s
          phytoplankton_mortality::FT = 5.8e-7, # 1/s
          small_detritus_remineralisation_rate::FT = 5.88e-7, # 1/s
          large_detritus_remineralisation_rate::FT = 5.88e-7, # 1/s
          phytoplankton_exudation_fraction::FT = 0.05,
          nitrification_rate::FT = 5.8e-7, # 1/s
          ammonia_fraction_of_exudate::FT = 0.75, 
          ammonia_fraction_of_excriment::FT = 0.5,
          ammonia_fraction_of_detritus::FT = 0.0,
          phytoplankton_redfield::FT = 6.56, # mol C/mol N
          organic_redfield::FT = 6.56, # mol C/mol N
          phytoplankton_chlorophyll_ratio::FT = 1.31, # g Chl/mol N
          organic_carbon_calcate_ratio::FT = 0.1, # mol CaCO₃/mol C
          respiration_oxygen_nitrogen_ratio::FT = 10.75, # mol O/molN
          nitrification_oxygen_nitrogen_ratio::FT = 2.0, # mol O/molN
          slow_sinking_mortality_fraction::FT = 0.5, 
          fast_sinking_mortality_fraction::FT = 0.5,
          dissolved_organic_breakdown_rate::FT = 3.86e-7, # 1/s
          zooplankton_calcite_dissolution::FT = 0.3,

          surface_photosynthetically_active_radiation = default_surface_PAR,

          light_attenuation_model::LA =
              TwoBandPhotosyntheticallyActiveRadiation(; grid, 
                                                         surface_PAR = surface_photosynthetically_active_radiation),
          sediment_model::S = nothing,

          carbonates::Bool = false,
          oxygen::Bool = false,
          variable_redfield::Bool = false,

          sinking_speeds = (sPOM = 3.47e-5, bPOM = 200/day),
          open_bottom::Bool = true,

          scale_negatives = false,

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

Construct an instance of the LOBSTER biogeochemical model.

Keyword Arguments

• grid: (required) the geometry to build the model on, required to calculate sinking
• phytoplankton_preference, ..., dissolved_organic_breakdown_rate: LOBSTER parameter values
• surface_photosynthetically_active_radiation: funciton (or array in the future) for the photosynthetically available radiation at the surface, should be shape f(x, y, t)
• light_attenuation_model: light attenuation model which integrated the attenuation of available light
• sediment_model: slot for AbstractSediment
• carbonates, oxygen, and variable_redfield: include models for carbonate chemistry and/or oxygen chemistry and/or variable redfield ratio dissolved and particulate organic matter
• 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=(3, 3, 30), extent=(10, 10, 200));

julia> model = LOBSTER(; grid)
LOBSTER{Float64} with carbonates ❌, oxygen ❌, variable Redfield ratio ❌ and (:sPOM, :bPOM) sinking 
 Light attenuation: Two-band light attenuation model (Float64)
 Sediment: Nothing
 Particles: Nothing
 Modifiers: Nothing

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³)
• Photosynthetically available radiation: PAR (einstein/m²/day)

Optional

• Ammonia: NH₄ (mmol N/m³)

SLatissima(; architecture :: AR = CPU(),
             growth_rate_adjustment :: FT = 4.5,
             photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60),
             minimum_carbon_reserve :: FT = 0.01,
             structural_carbon :: FT = 0.2,
             exudation :: FT = 0.5,
             erosion :: FT = 0.22,
             saturation_irradiance :: FT = 90 * day/ (10 ^ 6),
             structural_dry_weight_per_area :: FT = 0.5,
             structural_dry_to_wet_weight :: FT = 0.0785,
             carbon_reserve_per_carbon :: FT = 2.1213,
             nitrogen_reserve_per_nitrogen :: FT = 2.72,
             minimum_nitrogen_reserve :: FT = 0.0126,
             maximum_nitrogen_reserve :: FT = 0.0216,
             growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)),
             growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2,
             maximum_specific_growth_rate :: FT = 0.18,
             structural_nitrogen :: FT = 0.0146,
             photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24,
             photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24,
             photosynthesis_ref_temp_1 :: FT = 285.0,
             photosynthesis_ref_temp_2 :: FT = 288.0,
             photoperiod_1 :: FT = 0.85,
             photoperiod_2 :: FT = 0.3,
             respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24,
             respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24,
             respiration_ref_temp_1 :: FT = 285.0,
             respiration_ref_temp_2 :: FT = 290.0,
             photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1),
             photosynthesis_low_temp :: FT = 271.0,
             photosynthesis_high_temp :: FT = 296.0,
             photosynthesis_high_arrhenius_temp :: FT = 1414.87,
             photosynthesis_low_arrhenius_temp :: FT = 4547.89,
             respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1),
             current_speed_for_0p65_uptake :: FT = 0.03,
             nitrate_half_saturation :: FT = 4.0,
             ammonia_half_saturation :: FT = 1.3,
             maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6),
             maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6),
             current_1 :: FT = 0.72,
             current_2 :: FT = 0.28,
             current_3 :: FT = 0.045,
             respiration_reference_A :: FT = 1.11e-4 * 24,
             respiration_reference_B :: FT = 5.57e-5 * 24,
             exudation_redfield_ratio :: FT = Inf,

             prescribed_velocity :: U = 0.1,

             #position
             x :: P = on_architecture(architecture, [0.0])
             y :: P = on_architecture(architecture, zeros(Float64, length(x))),
             z :: P = on_architecture(architecture, zeros(Float64, length(x))),

             #properties
             A :: P = on_architecture(architecture, ones(Float64, length(x)) * 30),
             N :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.01),
             C :: P = on_architecture(architecture, ones(Float64, length(x)) * 0.1),

             #feedback
             nitrate_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))),
             ammonia_uptake :: P = on_architecture(architecture, zeros(Float64, length(x))),
             primary_production :: P = on_architecture(architecture, zeros(Float64, length(x))),
             frond_exudation :: P = on_architecture(architecture, zeros(Float64, length(x))),
             nitrogen_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))),
             carbon_erosion :: P = on_architecture(architecture, zeros(Float64, length(x))),

             custom_dynamics :: F = no_dynamics,

             scalefactor :: FT = 1.0,
             latitude :: FT = 57.5)

Keyword Arguments

• architecture: the architecture to adapt arrays to
• growth_rate_adjustment, ..., exudation_redfield_ratio: parameter values
• prescribed_velocity: functions for the relative velocity f(x, y, z, t)
• x,y and z: positions of the particles
• A, N, and C: area, nitrogen, and carbon reserves
• nitrate_uptake ... carbon_erosion: diagnostic values coupled to tracer fields
• custom_dynamics: place to add any function of form f!(particles, model, bgc, Δt)
• scalefactor: scalar scaling for tracer coupling
• latitude: model latitude for seasonal growth modulation

CarbonChemistryModel to solve chemical equilibrium parameterisations

CarbonChemistry(; ionic_strength = IonicStrength(),
                  solubility = K0(),
                  carbonic_acid = (K1 = K1(), K2 = K2()),
                  boric_acid = KB(),
                  water = KW(),
                  sulfate = KS(; ionic_strength),
                  fluoride = KF(; ionic_strength),
                  phosphoric_acid = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()),
                  silicic_acid = KSi(; ionic_strength),
                  calcite_solubility = KSP_calcite(),
                  density_function = teos10_polynomial_approximation)

Carbon chemistry model capable of solving for sea water pCO₂ from DIC and total alkalinity or DIC and pH.

Default form from Dickson, A.G., Sabine, C.L. and Christian, J.R. (2007), Guide to Best Practices for Ocean CO 2 Measurements. PICES Special Publication 3, 191 pp.

See each parameters documentation for origional sources.

Example

julia> using OceanBioME

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

julia> pCO₂ = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35)
1308.0843992121615

julia> pH = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35, return_pH = true)
7.502534641304366

julia> pCO₂_higher_pH = carbon_chemistry(; DIC = 2000, T = 10, S = 35, pH = 7.5)
1315.6558976217746

(p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing,
                       return_pH = false,
                       boron = 0.000232 / 10.811 * S / 1.80655,
                       sulfate = 0.14 / 96.06 * S / 1.80655,
                       fluoride = 0.000067 / 18.9984 * S / 1.80655,
                       silicate = 0,
                       phosphate = 0,
                       upper_pH_bound = 14,
                       lower_pH_bound = 0)

Calculates fCO₂ in sea water with DIC, 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.

Alternativly, pH is returned if return_pH is true.

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(; constant =  61.2172,
     inverse_T = -3633.86,
     log_T = -9.67770,
     S =  0.011555,
     S² = -0.0001152,
     pressure_correction = 
        PressureCorrection(; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877))

Parameterisation for aquious carbon dioxide - bicarbonate dissociation equilibrium constant.

CO₂*(aq) + H₂O ⇌ H₂CO₃ ⇌ HCO₃⁻ + H⁺

K₁ = [H⁺][HCO₃⁻]/[CO₂*]

Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119).

K2(; constant = -25.9290,
     inverse_T = -471.78,
     log_T =  3.16967,
     S =  0.01781,
     S² = -0.0001122,
     pressure_correction = 
        PressureCorrection(; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475))

Parameterisation for bicarbonate dissociation equilibrium constant.

HCO₃⁻ ⇌ CO₃²⁻ + H⁺

K₂ = [H⁺][CO₃²⁻]/[HCO₃⁻]

Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119).

KB(; constant =  148.0248,
     inverse_T = -8966.90,
     inverse_T_sqrt_S = -2890.53,
     inverse_T_S = -77.942,
     inverse_T_sqrt_S³ =  1.728,
     inverse_T_S² = -0.0996,
     sqrt_S = 137.1942,
     S = 1.62142,
     log_T = -24.4344,
     log_T_sqrt_S = -25.085,
     S_log_T = -0.2474,
     T_sqrt_S =  0.053105,  
     pressure_correction = 
        PressureCorrection(; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0))

Parameterisation for boric acid equilibrium with water.

B(OH)₃ + H₂O ⇌ B(OH)₄⁻ + H⁺

Kᵇ = [H⁺][B(OH)₄⁻]/[B(OH)₃]

Default values from Dickson (1990, Deep-Sea Res., 37, 755–766).

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(a₀, a₁, a₂,
                   b₀, b₁, b₂,
                   R)

Parameterisation for the pressure effect on thermodynamic constants.

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

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

alkalinity_residual(H, p)

Returns the difference between total alkalinity computed from H(hydrogen ion concentration),DIC,borate,sulfate,phosphate,silicate, andfluorideconcentration and chemical equilibrium constants specified inp, and the specified totalAlk`alinity.

TAlk = [HCO₃⁻] + 2[CO₃²⁻] + [B(OH)₄⁻] + [OH⁻] + [HPO₄²⁻] + 2[PO₄³⁻] + [SiO(OH)₃⁻] 
       + [NH₃] + [HS⁻] - [H⁺] - [HSO₄⁻] - [HF] - [H₃PO₄] + minor acids and bases

Concentrations diagnosed as specified in Dickson et. al best practice descried in CarbonChemistry docstring.

Note ammonia (NH₃) is not currently included.

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

TwoBandPhotosyntheticallyActiveRadiation(; grid, 
                                           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)

struct InstantRemineralisation

Hold the parameters and fields the simplest benthic boundary layer where organic carbon is assumed to remineralise instantly with some portion becoming N, and a fraction being permanently buried.

Burial efficiency from Dunne et al. (2007).

InstantRemineralisation(; grid,
    burial_efficiency_constant1::FT = 0.013,
    burial_efficiency_constant2::FT = 0.53,
    burial_efficiency_half_saturation::FT = 7)

Return a single-layer instant remineralisaiton model for NPZD bgc models.

Example

using OceanBioME, Oceananigans, OceanBioME.Sediments

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

sediment_model = InstantRemineralisation(; grid)

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

SimpleMultiG(; 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,
               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]))

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, OceanBioME.Sediments

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

julia> sediment_model = SimpleMultiG(; grid)
┌ Warning: Sediment models are an experimental feature and have not yet been validated.
└ @ OceanBioME.Models.Sediments ~/Documents/Projects/OceanBioME.jl/src/Models/Sediments/simple_multi_G.jl:104
[ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
Single-layer multi-G sediment model (Float64)

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(; carbon_chemistry = CarbonChemistry(),
                                            transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = CarbonDioxidePolynomialSchmidtNumber()),
                                            air_concentration = 413, # ppmv
                                            wind_speed = 2,
                                            water_concentration = nothing,
                                            silicate_and_phosphate_names = nothing,
                                            kwargs...)

Returns a FluxBoundaryCondition for the gas exchange between carbon dioxide dissolved in the water specified by the carbon_chemisty model, and air_concentration with transfer_velocity (see GasExchangeBoundaryCondition for details).

silicate_and_phosphate_names should either be nothing, a 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(; 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(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()),
                                     water_concentration = OxygenConcentration(),
                                     air_concentration = 9352.7, # mmolO₂/m³
                                     wind_speed = 2,
                                     kwagrs...)

Returns a FluxBoundaryCondition for the gas exchange between oxygen dissolved in the water specified by the the OxygenConcentration in the base model, and air_concentration with transfer_velocity (see GasExchangeBoundaryCondition for details).

kwargs are passed on to GasExchangeBoundaryCondition.

OxygenPolynomialSchmidtNumber(; 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(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(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(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(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(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(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(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(scale_factor = 0.27 / hour / 100)

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

Wanninkhof09(constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100)

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

Wanninkhof14(scale_factor = 0.251 / hour / 100)

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

Wanninkhof99(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 = (T = (t) -> 0, ))

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

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.