Biogeoceananigans.jl

[glwagner]

There's been a bunch of interest in using Oceananigans for physical-biogeochemical interaction studies --- problems where systems of reacting tracers that represent either oceanic biological systems, chemical reactions and cycles, or both interact with turbulence and fluid dynamics simulated by Oceananigans.

In order to organize the community's work (and also to support some work at Clima on biogeochemistry in Oceananigans), I propose that we create a new package that interfaces with Oceananigans --- could it be... Biogeoceananigans.jl... ? --- to:

1. Facilitate sharing code, and collaboration on implementation and testing of biogeochemistry models to be used in Oceananigans simulations
2. Develop documentation and a suite of examples to illustrate usage, setup, and analysis of numerical experiments with biogeochemistry

To achieve either of these there's no question we need a particular place to collaborate on re-useable code (rather than working independently). But also, I think the Oceananigans.jl repository is not the best repo to use to achieve the above goals, because it's big and complex, which might make it harder for potential developers to contribute and see their place. I also think it would slow development down, because, for example, we'll have to make sure all unit tests for differential operators pass before we can add a new biogeochemistry model implementation. I think development might be faster and more accessible if we do it in a different repo.

There are also some details to discuss regarding implementation. Oceananigans' design already supports reacting systems via Forcing. Oceananigans Forcing are arbitrary functions of spatial coordinates, time, prognostic model fields, and forcing function parameters --- or alternatively, indices i, j, k, grid, clock, and a NamedTuple of model fields that can be indexed into arbitrarily. @iuryt and @syou83syou83 (and perhaps others) have experimented in this direction.

However, I think we might benefit from adding a model property to both NonhydrostaticModel and HydrostaticFreeSurfaceModel that's specifically dedicated to biogeochemistry, and designing an interface that allows users / other packages (like Biogeoceananigans.jl) to build custom biogeochemistry types. One advantage I think is we'll have more flexibility in designing the code that hooks into Oceananigans than if we limit ourselves to using model.forcing. I also think it will yield a cleaner and more interpretable user API: biogeochemistry models are sometimes complicated, so it's helpful to disentangle their specification from the specification of more "vanilla" forcing terms, like sponge layers, etc. We may want special time-stepping methods for some types of biogeochemistry in the future.

I also think we can build an interface for biogeochemistry that uses similar or identical syntax to Forcing (like requiring functions of x, y, z, t, fields --- etc) because this will allow directly transferring code for biogeochemistry models that's first implemented as a system of forcing functions (a very common and useful development paradigm).

What do others think about this? @iuryt @johnryantaylor @rafferrari @simone-silvestri @christophernhill @jm-c @sandreza

[francispoulin]

This sounds like a big first step in the right direction. Having a another repo called biogeoceanigans.jl seems like a great choice! It would be nice it this was compatable with all the models.

[iuryt]

I draft a repo on this matter here.
I tried to summarize my ideas in the README:

Most of the interelationships between biogeochemical tracers could be implemented using the forcing functions (see Convecting Plankton example from Oceananigans documentation), but commonly the light or light-limiting growth of the phytoplankton are averaged over the mixed layer in order to parameterized the effects of the mixing in the upper boundary layer. In other words, the mixing timescale is usually longer than photosynthesis, but shorter than cell-division time, which makes phytoplankton see the average light or grow at the average light-limiting growth rate while on the mixed layer. Other example of processes that are not easily implemented on Oceananigans include the shading of phytoplankton to deeper layers and the sinking velocity. The latter is currently being implemented on Oceananigans by this Pull Request with a new forcing called AdvectiveForcing.

For now, this repository basically gives a set of modules that can be used to estimate the mixed-layer depth, estimate the phytoplankton shading and calculate the light-limiting growth, which is defined by the user as a function, and then the light (or the light-limiting growth) can be averaged over the mixed layer.

So basically the repository can simply give all auxiliary functions that are needed for the biogeochemical models, but not the interelationship functions, which can be defined as forcing functions and added to a list of examples on the documentation.

I don't mind giving up this repository and merging or adding it to the organization if this is somehow a good starting point.
I am also interested in seeing what the others think about that.

[glwagner]

@iuryt that's a beautiful summary! That repository is a great starting point! All we need is the word geo I think... 😄

Reading your description I'm realizing there's another point to be made about "embedded" biogeochemistry versus "external" when the biogeochemistry has it's own diagnostic state that must be precomputed before a time-step (the example of light growth is a good one). Embedded helps with this, because users don't have to add callbacks, etc themselves.

[iuryt]

Yeah.. I am not yet sure how this would work embedded. By now, what I have is just some functions that the user has to add as callbacks.
e.g.
https://github.com/iuryt/Bioceananigans.jl/blob/ab35a7235fea532fe7ed283ce6e06ac12f378dae/test/NNP.jl#L109-L117
btw.. this test script is not yet fully developed 😄.

Thinking more, it would be so the function automatically adds the callbacks to the simulation?

const Δb=(g/ρₒ) * 0.03
h = Field{Center, Center, Nothing}(grid) 
MixedLayerDepth(simulation, h, Δb)


# constants for the NP model
const μ₀ = 1/day   # surface growth rate
const m = 0.015/day # mortality rate due to virus and zooplankton grazing
const Kw = 0.059 # meter^-1
const kn = 0.75
const kr = 0.5

#  https://doi.org/10.1029/2017GB005850
const chl2c = 0.06 # average value for winter in North Atlantic

const α = 0.0538/day

const average = :growth
const shading = true

light_growth = Field{Center, Center, Center}(grid)

# time evolution of shortwave radiation (North Atlantic)
@inline Lₒ(t) = 116.5 * sin( 2π * ( t / days + 50 ) / 375.7 + 1.3 ) + 132.3
# evolution of the available light at the surface
@inline light_function(t, z) = 0.43 * Lₒ(t) * exp( z * Kw )
# light profile
@inline light_growth_function(light) = μ₀ * ( light * α ) / sqrt( μ₀^2 + ( light * α )^2 )

phytoplankton = ("P",) # list of phytoplankton tracers  
LightGrowth(light_growth, h, simulation, phytoplankton, light_function, light_growth_function, time(simulation), average, shading, chl2c)

[christophernhill]

@iuryt https://github.com/iuryt/Bioceananigans.jl looks like an awesome start!!! 🚀

@glwagner arguably it could be biogeochemonanigans.jl or biogeocheminanigans.jl ?

Would it be interesting to think about an interface that could be helpful for projects like https://github.com/openjournals/joss-reviews/issues/4207/https://github.com/JuliaOcean/PlanktonIndividuals.jl from @zhenwu0728 and @gaelforget too at some point?

Note -

1. in the end both Eulerian and Lagrangian formulations of biogeocheminanigans/biogeochemonanigans/Bioceananigans/Bichemoceananigans style stuff end up as mostly pointwise functions. So a number of bits of an API framework could be common across Eulerian and Lagrangian I think.

2. On the fluid solver side the API could also support reactive and passive transport in an atmosphere (i.e. work with atmoscheminanigans 😄 )

[glwagner]

@glwagner arguably it could be biogeochemonanigans.jl or biogeocheminanigans.jl ?

I know chemistry! Ahaha I like those. Biogeochemmingaround.jl? Ok, let's brainstorm.

Thinking more, it would be so the function automatically adds the callbacks to the simulation?

There are multiple possible places where one might update the biogeochemical_state. An "external" paradigm requires Callback. But an "embedded" paradigm can hard-code functionality into update_state(model)!:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl

I think we want to pursue code design that allows these paradigms to be interchanged. I think it's crucial that functionality developed "externally" can be "copy/pasted" into Biogeoceananigans.jl. This probably would be less pretty than it sounds (style, code quality, names...), but I think is still a good organizing principle for the design.

Would it be interesting to think about an interface that could be helpful for projects like https://github.com/openjournals/joss-reviews/issues/4207/https://github.com/JuliaOcean/PlanktonIndividuals.jl from @zhenwu0728 and @gaelforget too at some point?

I see very much what you're saying. There's an opportunity to also support biogeochemistry on Lagrangian particles.

[iuryt]

That's an awesome idea @christophernhill !

I haven't yet used a Lagrangian formulation for biogeochemical models, but this it would be cool to have some way to support both on this.

The name Oceananigans is already a challenge for me to pronounce as a non-native speaker. I must confess I had to practice in the front of the mirror many times. haha

What about Biogeochemigans.jl?

We can list some of the options and make a poll. Is it possible to create a poll on github?

[rafferrari]

I suggest we wait on this for the time being. We are in discussions to start a formal collaboration to develop the ocean biogeochemistry component of the CliMA model. We need to plan our strategy with the new collaborators who will join us. It will happen soon and we can then have a discussion with all people involved on what is the best strategy moving forward.

[francispoulin]

That sounds like a great idea @rafferrari . I would be happy to participate in this in the future if you are looking for people to get involved.

[gaelforget]

The more the merrier! Contributors are generally welcome on our end too. Links @christophernhill in his comment were broken. So just in case:

• JOSS issue (accepted), paper (proofs)
• Github repo (JuliaOcean)

PlanktonIndividuals.jl: A GPU supported individual-based phytoplankton life cycle model

[glwagner]

The more the merrier! Contributors are generally welcome on our end too. Links @christophernhill in his comment were broken. So just in case:

• JOSS issue (accepted), paper (proofs)
• Github repo (JuliaOcean)

PlanktonIndividuals.jl: A GPU supported individual-based phytoplankton life cycle model

@gaelforget it looks like you might have some Oceananigans examples which is super cool!! Can you point to them for us? Have you considered expanding to "online" capabilities? What are the pros and cons of that?

[glwagner]

I suggest we wait on this for the time being. We are in discussions to start a formal collaboration to develop the ocean biogeochemistry component of the CliMA model. We need to plan our strategy with the new collaborators who will join us. It will happen soon and we can then have a discussion with all people involved on what is the best strategy moving forward.

I agree, this issue is for planning purposes, but there's no timeline for setting up a new package. I think the package itself is easy, but designing the interface that allows different "biogeochemical models" to plug in to Oceananigans will require a bit more care. On the Oceananigans side, it seems like some necessary features of a biogeochemistry interface are

• User-specified biogeochemistry model "automatically" adds the tracers it needs in the model constructor (does not require user to explicitly specify tracers): model_tracers = add_biogeochemical_tracers(biogeochemistry, user_tracers).
• There's a "biogeochemical tendency" term in the tracer equations that's controlled by model.biogeochemistry type.
• There's an interface in update_state!(model) like update_biogeochemical_state!(model.biogeochemistry, model) so that biogeochemistry models can precompute a state (eg cumulative integrals to determine light penetration) when needed.

Perhaps further consultation will produce additional requirements. One big one that I see is for chemistry models that require subcycling or special numerical methods for time-stepping (ie stiff carbonate chemistry systems as in Smith et al 2018). Supporting those cases will be more challenging.

I think our Lagrangian particles already (or are intended to) support dynamics-on-particles so there may not be much to do there on the Oceananigans side.

[gaelforget]

The more the merrier! Contributors are generally welcome on our end too. Links @christophernhill in his comment were broken. So just in case:

• JOSS issue (accepted), paper (proofs)
• Github repo (JuliaOcean)

PlanktonIndividuals.jl: A GPU supported individual-based phytoplankton life cycle model

@gaelforget it looks like you might have some Oceananigans examples which is super cool!! Can you point to them for us? Have you considered expanding to "online" capabilities? What are the pros and cons of that?

The two models we currently have some support for (offline) are MITgcm and oceananigans.jl.

I have examples for both models in https://gaelforget.github.io/ClimateModels.jl/dev/examples/ which you can run the JuliaClimate sandbox https://juliaclimate.github.io/Notebooks/ and generate the kind of model output that PlanktonIndividuals.jl or IndividualDisplacement.jl can then ingest as input for offline computation. The PlanktonIndividuals.jl documentation has corresponding examples on the ecology + bgc side.

An online mode with models like MITgcm and oceananigans.jl would be useful in my view, in the case of both PlanktonIndividuals.jl or IndividualDisplacement.jl. With this being said, I tend to most like the idea of small'ish modeling components developed outside of monolitic modeling efforts that emphasize code re-usability and integration with multiple models.

For completeness, here is the reference to the JOSS paper on IndividualDisplacement.jl -- https://joss.theoj.org/papers/10.21105/joss.02813.pdf -- and the GitHub repo

[glwagner]

I guess if you have a (julia) function that takes a particle time-step, given a velocity field (ie 3D array), then you can use PlanktonIndividuals.jl online with Oceananigans now, using the Callback Oceananigans feature that's designed for integration with other models. So maybe it already works in fact.

[iuryt]

For now, I will change the name of the package to Biogeochemigans.jl and wait till further discussions have advanced on this matter. I will be using this on a paper under review, so I have to move forward anyway. As I said, I don't mind changing my focus to another repository or using part of the work I did for this one (if it is somehow useful). Thanks for all this discussion, I am glad to help with this project.

[glwagner]

We're close to having this, I'm going to move this to a discussion.