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)

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.

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,

                                           particles::P = 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
• particles: slot for BiogeochemicalParticles

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::SPAR = 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 = false,

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

          particles::P = 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
• particles: slot for BiogeochemicalParticles

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 = arch_array(architecture, [0.0])
             y :: P = arch_array(architecture, zeros(Float64, length(x))),
             z :: P = arch_array(architecture, zeros(Float64, length(x))),

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

             #feedback
             nitrate_uptake :: P = arch_array(architecture, zeros(Float64, length(x))),
             ammonia_uptake :: P = arch_array(architecture, zeros(Float64, length(x))),
             primary_production :: P = arch_array(architecture, zeros(Float64, length(x))),
             frond_exudation :: P = arch_array(architecture, zeros(Float64, length(x))),
             nitrogen_erosion :: P = arch_array(architecture, zeros(Float64, length(x))),
             carbon_erosion :: P = arch_array(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

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

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 = (x, y, t) -> 100 * max(0.0, cos(t * π / 12hours)))

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)

Boundary conditions for air/sea and sediment flux.

Currently implemented:

• gasexchange (Wanninkhof, 1992)

  • Generic air-sea flux model described by Wanninkhof (1992) but only setup for CO₂ and O₂.
  • Forces the DIC and oxygen fields, and requires temp (in centigrade) and salinity, plus current DIC and ALK concentration.

• Sediments

  • Soetaert (Soetaert et al., 2000)
    • simple (integrated) sediment model described by Soetaert et al. (2000), where organic matter that sinks to the bottom is stored, decays into NO₃ and NH₄, and takes up O₂ in the process.
    • Extended to attribute the corresponding release of DIC.
    • Forced by O₂, NO₃, NH₄ and particle concentration in bottom cell.
  • Instant remineralisation
    • simple model from Aumont et al. (2015), where sinking organic matter is instantly remineralised and returned to the bottom cell
    • some fraction is stored permanently in the sediment at an efficiency given by Dunne et al. (2007)

GasExchange(; gas,
              schmidt_params::ScP = (CO₂ = (A=2073.1, B=125.62, C=3.6276, D=0.043219),
                               O₂ = (A=1953.4, B=128.0, C=3.9918, D=0.050091))[gas],
              solubility_params::βP = (CO₂ = (A₁=-60.2409, A₂=93.4517, A₃=23.3585, B₁=0.023517, B₂=-0.023656, B₃=0.0047036),
                                 O₂ = (A₁=-58.3877, A₂=85.8079, A₃=23.8439, B₁=-0.034892, B₂=0.015568, B₃=-0.0019387))[gas],
              ocean_density::FT = 1026, # kg/m³
              air_concentration::AC = (CO₂ = 413.4, O₂ = 9352.7)[gas], # ppmv, mmolO₂/m³ (20.95 mol O₂/mol air, 0.0224m^3/mol air)
              air_pressure::FT = 1.0, # atm
              average_wind_speed::FT = 10, # m/s
              field_dependencies = (CO₂ = (:DIC, :ALK), O₂ = (:OXY, ))[gas])

Construct an Oceananigans FluxBoundaryCondition for the exchange of gas with the relevant tracer (i.e., DIC for CO₂ and oxygen for O₂). Please see note for other gases.

Keyword arguments

• gas: (required) the gas to be exchanged, if :CO₂ or :O₂ are specified then all other settings may be infered
• schmidt_params : named tuple of parameters for calculating the Schmidt number using the parameterisation of Wanninkhof (1992)
• solubility_params : named tuple of parameters for calculating the solubility (for O₂ the Bunsen solubility and CO₂ K₀, see note)
• ocean_density : density of the ocean in kg/m³
• air_concentratio : concentration of the gas in air in relivant units (i.e. ppmv for CO₂ and mmol O₂/m³ for O₂), can also be a function of x, y, t, or a field
• air_pressure : air pressure in atm (only used for CO₂), can also be a function of x, y, t, or a field
• average_wind_speed : average wind speed at 10m used to calculate the gas transfer velocity by the Wanninkhof (1992) parameterisation
• field_dependencies : tracer fields that gas exchange depends on, if the defaults have different names in your model you can specify as long as they are in the same order
• pCO₂ : pCO₂ calculator

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 = arch_array(architecture(grid), [- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232]),
               denitrification_params = arch_array(architecture(grid), [- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662]),
               anoxic_params = arch_array(architecture(grid), [- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143]),
               solid_dep_params = arch_array(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.Boundaries.Sediments ~/OceanBioME.jl/src/Boundaries/Sediments/simple_multi_G.jl:102
[ Info: This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC.
Single-layer multi-G sediment model (Float64)

Integrate biogeochemical models on a single point

BoxModel(; biogeochemistry::B,
           stop_time::FT = 0.0,
           forcing = NamedTuple(),
           timestepper::TS = RungeKutta3TimeStepper((required_biogeochemical_tracers(biogeochemistry)..., required_biogeochemical_auxiliary_fields(biogeochemistry)...)),
           Δt::FT = 1.0,
           clock::C = Clock(0.0, 0, 1))

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

Keyword Arguments

• biogeochemistry: (required) an OceanBioME biogeochemical model, most models must be passed a grid which can be set to BoxModelGrid() for box models
• stop_time: end time of simulation
• forcing: NamedTuple of additional forcing functions for the biogeochemical tracers to be integrated
• timestepper: Timestepper to integrate model, only available is currently RungeKutta3TimeStepper
• Δt: time step length
• clock: Oceananigans clock to keep track of time

SaveBoxModel(filepath::FP)

Construct object to save box model outputs at filepath.

Arguments

• filepath - path to save results to

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

run!(model::BoxModel; feedback_interval = 1000, save_interval = Inf, save = nothing)

Run a box model.

Arguments

• model - the BoxModel to solve

Keyword Arguments

• feedback_interval: how often (number of iterations) to display progress
• save_interval: how often (number of iterations) to save output
• save: SaveBoxModel object to specify how to save output

TODO: should abstract out to simulation like Oceananigans to add e.g. callbacks