# Implementing a new models

Here we describe how OceanBioME defines biogeochemical (BGC) models, how this varies from Oceananigans, and how to implement your own model.

## Model structure

OceanBioME BGC models are `struct`

s of type `ContinuousFormBiogeochemistry`

, which is of abstract type `AbstractContinuousFormBiogeochemistry`

from Oceananigans. In Oceananigans this describes BGC models which are defined using continuous functions (depending continuously on $x$, $y$, and $z$) rather than discrete functions (depending on $i$, $j$, $k$). This allows the user to implement the BGC model equations without worrying about details of the grid or discretization, and then Oceananigans handles the rest.

OceanBioME's `ContinuousFormBiogeochemistry`

adds a layer on top of this which makes it easy to add light attenuation models, sediment, and biologically active particles (or individual-based models). OceanBioME's `ContinuousFormBiogeochemistry`

includes parameters in which the types of these components are stored. This means that these model components will automatically be integrated into the BGC model without having to add new methods to call Oceananigans functions.

## Implementing a model

The nature of multiple dispatch in Julia means that we define new BGC models as new types. You can then define methods to this type which are used by OceanBioME and Oceananigans to integrate the model.

### The basics

For this example we are going to implement the simple Nutrient-Phytoplankton model similar to that used in (Chen *et al.*, 2015), although we neglect the nutrient in/outflow terms since they may be added as boundary conditions, and modified to conserve nitrogen.

The first step is to import the abstract type from OceanBioME, some units from Oceananigans (for ease of parameter definition), and `import`

some functions from Oceananigans in order to add methods to:

```
using OceanBioME: Biogeochemistry
using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry
using Oceananigans.Units
import Oceananigans.Biogeochemistry: required_biogeochemical_tracers,
required_biogeochemical_auxiliary_fields,
biogeochemical_drift_velocity
```

We then define our `struct`

with the model parameters, as well as slots for the particles, light attenuation, and sediment models:

```
@kwdef struct NutrientPhytoplankton{FT, W} <: AbstractContinuousFormBiogeochemistry
base_growth_rate :: FT = 1.27 / day # 1 / seconds
nutrient_half_saturation :: FT = 0.025 * 1000 / 14 # mmol N / m³
light_half_saturation :: FT = 300.0 # micro einstein / m² / s
temperature_exponent :: FT = 0.24 # 1
temperature_coefficient :: FT = 1.57 # 1
optimal_temperature :: FT = 28.0 # °C
mortality_rate :: FT = 0.15 / day # 1 / seconds
crowding_mortality_rate :: FT = 0.004 / day / 1000 * 14 # 1 / seconds / mmol N / m³
sinking_velocity :: W = 2 / day
end
```

Here, we use descriptive names for the parameters. Below, each of these parameters correspond to a symbol (or letter) which is more convenient mathematically and when defining the BGC model functions. In the above code we used `@kwdef`

to set default values for the models so that we don't have to set all of these parameters each time we use the model. The default parameter values can optionally be over-ridden by the user when running the model. We have also included a `sinking_velocity`

field in the parameter set to demonstrate how we can get tracers (e.g. detritus) to sink. We also need to define some functions so that OceanBioME and Oceananigans know what tracers and auxiliary fields (e.g. light intensity) we use:

```
required_biogeochemical_tracers(::NutrientPhytoplankton) = (:N, :P, :T)
required_biogeochemical_auxiliary_fields(::NutrientPhytoplankton) = (:PAR, )
```

`required_biogeochemical_auxiliary_fields (generic function with 5 methods)`

Next, we define the functions that specify how the phytoplankton $P$ evolve. In the absence of advection and diffusion (both of which are handled by Oceananigans), we want the phytoplankton to evolve at the rate given by:

\[\frac{\partial P}{\partial t} = \mu g(T) f(N) h(PAR) P - mP - bP^2,\]

where $\mu$ corresponds to the parameter `base_growth_rate`

, $m$ corresponds to the parameter `mortality_rate`

, and $b$ corresponds to the parameter `crowding_mortality_rate`

. Here, the functions $g$, $f$, and $h$ are defined by:

\[\begin{align} g(T) &= c_1\exp\left(-c_2|T - T_{opt}|\right),\\ f(N) &= \frac{N}{k_N + N},\\ h(PAR) &= \frac{PAR}{k_P + PAR}, \end{align}\]

where $c_1$ corresponds to `temperature_coefficient`

, $c_2$ corresponds to `temperature_exponent`

, $T_{opt}$ corresponds to `optimal_temperature`

, $k_N$ corresponds to `nutrient_half_saturation`

, and $k_P$ corresponds to `light_half_saturation`

.

We turn this into a function for our model by writing:

```
@inline function (bgc::NutrientPhytoplankton)(::Val{:P}, x, y, z, t, N, P, T, PAR)
μ = bgc.base_growth_rate
m = bgc.mortality_rate
b = bgc.crowding_mortality_rate
growth = μ * g(bgc, T) * f(bgc, N) * h(bgc, PAR) * P
death = m * P + b * P ^ 2
return growth - death
end
@inline function g(bgc, T)
c₁ = bgc.temperature_coefficient
c₂ = bgc.temperature_exponent
Tₒ = bgc.optimal_temperature
return c₁ * exp(-c₂ * abs(T - Tₒ))
end
@inline function f(bgc, N)
kₙ = bgc.nutrient_half_saturation
return N / (N + kₙ)
end
@inline function h(bgc, PAR)
kₚ = bgc.light_half_saturation
return PAR / (PAR + kₚ)
end
```

`h (generic function with 1 method)`

The first parameter `::Val{:P}`

is a special value type that allows this function to be dispatched when it is given the value `Val(:P)`

. This is how Oceananigans tells the model which forcing function to use. At the start of the `NutrientPhytoplankton`

function we unpack some parameters from the model, then calculate each term, and return the total change (the gain minus the loss).

For this model, the nutrient evolution can be inferred from the rate of change of phytoplankton. Since this is a simple two variable model and the total concentration is conserved,

\[\frac{\partial N}{\partial t} = - \frac{\partial P}{\partial t}\]

Hence, we define the nutrient forcing using as the negative of the phytoplankton forcing

`@inline (bgc::NutrientPhytoplankton)(::Val{:N}, args...) = -bgc(Val(:P), args...)`

Finally, we need to define some functions to allow us to update the time-dependent parameters (in this case the PAR and temperature, $T$):

```
using OceanBioME: BoxModel
import OceanBioME.BoxModels: update_boxmodel_state!
function update_boxmodel_state!(model::BoxModel{<:NutrientPhytoplankton, <:Any, <:Any, <:Any, <:Any, <:Any})
getproperty(model.values, :PAR) .= model.forcing.PAR(model.clock.time)
getproperty(model.values, :T) .= model.forcing.T(model.clock.time)
end
```

`update_boxmodel_state! (generic function with 4 methods)`

Now we can run an example similar to the LOBSTER box model example:

```
using OceanBioME, Oceananigans.Units
const year = years = 365days
@inline PAR⁰(t) = 500 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
z = -10 # specify the nominal depth of the box for the PAR profile
@inline PAR(t) = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay
@inline temp(t) = 2.4 * cos(t * 2π / year + 50days) + 26
biogeochemistry = NutrientPhytoplankton()
model = BoxModel(; biogeochemistry, forcing = (; PAR, T = temp))
model.Δt = 5minutes
model.stop_time = 5years
set!(model, N = 15, P = 15)
# ## Run the model (should only take a few seconds)
@info "Running the model..."
run!(model, save_interval = 100, feedback_interval = Inf, save = SaveBoxModel("box_np.jld2"))
```

```
[ Info: Running the model...
[ Info: Reached 0 seconds
```

We can then visualise this:

```
using JLD2
vars = (:N, :P, :T, :PAR)
file = jldopen("box_np.jld2")
times = parse.(Float64, keys(file["values"]))
timeseries = NamedTuple{vars}(ntuple(t -> zeros(length(times)), length(vars)))
for (idx, time) in enumerate(times)
values = file["values/$time"]
for tracer in vars
getproperty(timeseries, tracer)[idx] = values[tracer]
end
end
close(file)
# ## And plot
using CairoMakie
fig = Figure(resolution = (1200, 480), fontsize = 20)
axN= Axis(fig[1, 1], ylabel = "Nutrient \n(mmol N / m³)")
lines!(axN, times / year, timeseries.N, linewidth = 3)
axP = Axis(fig[1, 2], ylabel = "Phytoplankton \n(mmol N / m³)")
lines!(axP, times / year, timeseries.P, linewidth = 3)
axPAR= Axis(fig[2, 1], ylabel = "PAR (einstein / m² / s)", xlabel = "Time (years)")
lines!(axPAR, times / year, timeseries.PAR, linewidth = 3)
axT = Axis(fig[2, 2], ylabel = "Temperature (°C)", xlabel = "Time (years)")
lines!(axT, times / year, timeseries.T, linewidth = 3)
fig
```

So now we know it works.

### Phytoplankton sinking

Now that we have a fully working BGC model we might want to add some more features. Another aspect that is easy to add is negative buoyancy (sinking). To-do this all we do is add a method to the Oceananigans function `biogeochemical_drift_velocity`

, and we use `::Val{:P}`

to specify that only phytoplankton will sink. Above, we set the default value of the parameter `bgc.sinking_velocity`

. We can override this when we call the BGC model like `NutrientPhytoplankton(; light_attenuation_model, sinking_velocity = 1/day)`

. Note that before using `biogeochemical_drift_velocity`

, we need to import several `Fields`

from Oceananigans:

```
using Oceananigans.Fields: ZeroField, ConstantField
biogeochemical_drift_velocity(bgc::NutrientPhytoplankton, ::Val{:P}) =
(u = ZeroField(), v = ZeroField(), w = bgc.sinking_velocity)
```

`biogeochemical_drift_velocity (generic function with 8 methods)`

### Sediment model coupling

Another aspect that OceanBioME includes is sediment models. Doing this varies between sediment models, but for the most generic and simplest, all we need to do is add methods to two functions:

```
using OceanBioME.Boundaries.Sediments: sinking_flux
import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers
@inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) =
sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers)
@inline carbon_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56
@inline remineralisation_receiver(::NutrientPhytoplankton) = :N
@inline sinking_tracers(::NutrientPhytoplankton) = (:P, )
```

`sinking_tracers (generic function with 5 methods)`

### Putting it together

Now that we have added these elements we can put it together into another simple example:

```
using Oceananigans, OceanBioME
using OceanBioME.Sediments: InstantRemineralisation
# define some simple forcing
@inline surface_PAR(t) = 200 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
@inline ∂ₜT(z, t) = - 2π / year * sin(t * 2π / year + 50days)
@inline κₜ(z, t) = 1e-2 * (1 + tanh((z - 50) / 10)) / 2 + 1e-4
# define the grid
grid = RectilinearGrid(topology = (Flat, Flat, Bounded), size = (32, ), extent = (100, ))
# setup the biogeochemical model
light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR)
sediment = InstantRemineralisation(; grid)
sinking_velocity = ZFaceField(grid)
w_sink(z) = 2 / day * tanh(z / 5)
set!(sinking_velocity, w_sink)
negative_tracer_scaling = ScaleNegativeTracers((:N, :P))
biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity);
light_attenuation,
sediment,
modifiers = negative_tracer_scaling)
# put the model together
model = NonhydrostaticModel(; grid,
biogeochemistry,
closure = ScalarDiffusivity(ν = κₜ, κ = κₜ),
forcing = (; T = ∂ₜT))
set!(model, P = 0.01, N = 15, T = 28)
# run
simulation = Simulation(model, Δt = 9minutes, stop_time = 1years)
simulation.output_writers[:tracers] = JLD2OutputWriter(model, model.tracers,
filename = "column_np.jld2",
schedule = TimeInterval(1day),
overwrite_existing = true)
simulation.output_writers[:sediment] = JLD2OutputWriter(model, model.biogeochemistry.sediment.fields,
indices = (:, :, 1),
filename = "column_np_sediment.jld2",
schedule = TimeInterval(1day),
overwrite_existing = true)
run!(simulation)
```

```
┌ Warning: Sediment models are an experimental feature and have not yet been validated
└ @ OceanBioME.Boundaries.Sediments ~/builds/kelp-6/oceanbiome/oceanbiome/src/Boundaries/Sediments/instant_remineralization.jl:57
[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (675.997 ms)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.936 seconds).
[ Info: Simulation is stopping after running for 56.097 seconds.
[ Info: Simulation time 365 days equals or exceeds stop time 365 days.
```

We can then visualise this:

```
N = FieldTimeSeries("column_np.jld2", "N")
P = FieldTimeSeries("column_np.jld2", "P")
sed = FieldTimeSeries("column_np_sediment.jld2", "N_storage")
fig = Figure()
axN = Axis(fig[1, 1], ylabel = "z (m)")
axP = Axis(fig[2, 1], ylabel = "z (m)")
axSed = Axis(fig[3, 1:2], ylabel = "Sediment (mmol N / m²)", xlabel = "Time (years)")
_, _, zc = nodes(grid, Center(), Center(), Center())
times = N.times
hmN = heatmap!(axN, times ./ year, zc, N[1, 1, 1:grid.Nz, 1:end]',
interpolate = true, colormap = Reverse(:batlow))
hmP = heatmap!(axP, times ./ year, zc, P[1, 1, 1:grid.Nz, 1:end]',
interpolate = true, colormap = Reverse(:batlow))
lines!(axSed, times ./ year, sed[1, 1, 1, :])
Colorbar(fig[1, 2], hmN, label = "Nutrient (mmol N / m³)")
Colorbar(fig[2, 2], hmP, label = "Phytoplankton (mmol N / m³)")
fig
```

We can see in this that some phytoplankton sink to the bottom, and are both remineralized back into nutrients and stored in the sediment.

### Final notes

When implementing a new model we recommend following a testing process as we have here, starting with a box model, then a column, and finally using it in a realistic physics scenarios. We have found this very helpful for spotting bugs that were proving difficult to decipher in other situations. You can also add `Individuals`

, light attenuation models, and sediment models in a similar fashion.