Add single pool soil carbon/biogeochemistry process

[BernhardAhrens]

Work in progress! Quick one pool model with a Q10 temperature sensitivity as an exercise for me to learn Terrarium. So far, I have to say I really like the package and was easy to get into with Max's example.

I was struggleing with the units of the solver. Is it unitless and just depends on the units of the rate constants (yr^-1 or s^-1)?

I also did not figure out the vertical dimension and how to implement advective and diffusive fluxes.

I wil hopefully come back to this later...

[bgroenks96]

Hey @BernhardAhrens! Thanks, this looks like a great start! I can help with some of these things. Do you mind if I also push to this branch?

[bgroenks96]

Also if you have a reference or some LaTeX lying around detailing the equation you want to implement, that would be super helpful!

[BernhardAhrens]

Yes, please feel free to push to the branch @bgroenks96

[BernhardAhrens]

Here is the supplement to my paper which contains all the equations: https://nextcloud.bgc-jena.mpg.de/s/ykwNWsErTnATtD3 There are a lot of pools you do not need to worry about. I would be interested how to implement the advection and diffusive fluxes (bioturbation and possibly cryoturbation, and possibly a prime candidate for hybrid modelling).

Here is the link to the whole paper: https://www.sciencedirect.com/science/article/abs/pii/S0038071720302091
Most microbial-mineral models converge to a structure similar to our model, so it would not be a problem to implement other models.

[bgroenks96]

Hi @BernhardAhrens, sorry for the long delay. I was busy working on some paper revisions.

Thanks for the links to your papers. These equations all look very doable, at least at first glance.

I would start with the diffusion terms since this is probably the easiest. You can take a look at the soil energy balance implementation to see how this is done with Oceananigans operators:

https://github.com/TUM-PIK-ESM/Terrarium.jl/blob/main/src/processes/soil/energy/soil_energy.jl#L118

Let me know if you have any questions.

Regarding the advection velocity $\omega(z)$, how is this actually calculated? Is it a prognostic variable with its own tendency? or is it constant in time and only variable in depth?

[BernhardAhrens]

There are two advection velocities, one advection velocity is just the water flux between layers. This advection velocity could come from your soil physics model - it was water fluxes from JSBACH site level or global runs as forcing in the past.

The $\omega(z)$ is similar to sedimentation velocity in marine models which pushes the previous profile [t-1] downwards with the thickness (m) of the new sediment per dt (yr-1). With this you end up with an advection velocity. In soils, we use the litterfall (kg m-2) with bulk density of organic matter (kg m-3) to calculate an advection velocity as well. Compared to sediments, we have two two extra $\omega$ terms in soils: one for carbon inputs via roots and one for the decomposition loss.

So, I think it is not a prognostic variable but derived from other fluxes. It varies with time and depth.

You can find more in this Appendix with a little sketch. In the end it might be a bit tricky to implement as we have to sum all inputs via roots and losses of SOC up to a depth $z$. Not a top priority 😉. This representation is relevant for peatlands and other organic carbon rich soils.

I'll have a look at the soil energy balance model and see where I get stuck.

[bgroenks96]

Hi @BernhardAhrens, sorry for the long silence. I hope your year is off to a good start!

Have you made any progress with this? After I finish some other tasks on getting a v0.1 done, I would probably circle back to this some time next month. I could prepare a skeleton implementation including also the of the Oceananigans operators and then you could take it from there?

[BernhardAhrens]

Hi @bgroenks96 Brian, I did not find time so far to work on it further but I plan in it soon and should have more time the next weeks. It would be really helpful if you showed the Oceananigans operators for transport. Lazaro @lazarusA was also looking at Oceananigans last week so could be starting point for us.

[bgroenks96]

Hi again @BernhardAhrens. I got a bit carried away in e4ecb34 and wrote more than just skeleton code, hope that's ok 😅. I adapted your initial implementation to follow the general style of process implementations in Terrarium and also to use kernels for GPU compatibility. I also moved your top-level code to a script in examples/scripts.

Note that the method interface for processes changed since we introduced it to you in November (see #68). Most notably, we no longer pass state into kernels because we found this to be a performance bottleneck. Instead, you have to explicitly collect the input and output fields in compute_auxiliary! and pass them to the kernel. We also made processes independent of models, so they have to accept as arguments any additional processes which they depend on for coupling.

This is very much still WIP of course, but I hope this gives you a sense of how Terrarium/Oceananigans works. Feel free to leave comments on the changes in this PR if you have any questions.

We are also working on adding a lot more documentation and we'll let you know once that's live so you can take a look.

Edit: Also, @BernhardAhrens, just to clarify, does "one pool" here mean one soil layer or one particular source/sink of carbon? like wood, microbial, etc. I assume it's the latter, but I just want to make sure I did the right thing by making the density_soc prognostic variable discretized over soil layers (XYZ).

[bgroenks96]

One other physical detail that has become apparent to me now is that we will definitely need to create state variables for soil component densities. The current approach (borrowed from CryoGrid) of pre-scaling the soil properties (like heat capacity) by constant densities doesn't make sense when you have prognostic carbon density because of course then the overall density becomes variable as well. I suppose we will also need to account for water advection due to the change in porosity? e.g. if carbon density decreases, and thus the bulk porosity changes, the total water mass in the soil volume would also change and that water needs to be routed to/from somewhere, otherwise we violate the water mass balance.

[bgroenks96]

Hi @BernhardAhrens, just wanted to check in and see if you had time to look at the initial implementation. Cheers!