Negative Tracers

[jagoosw]

I was looking through NEMO's PISCES source code and noticed that they enforce that tracers can't be negative with this routine:

         ! Handling of the negative concentrations
         ! The biological SMS may generate negative concentrations
         ! Trends are tested at each grid cell. If a negative concentrations 
         ! is created at a grid cell, all the sources and sinks at that grid 
         ! cell are scale to avoid that negative concentration. This approach 
         ! is quite simplistic but it conserves mass.

          xnegtr(:,:,:) = 1.e0
         DO jn = jp_pcs0, jp_pcs1
            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk)
               IF( ( tr(ji,jj,jk,jn,Kbb) + tr(ji,jj,jk,jn,Krhs) ) < 0.e0 ) THEN
                  ztra             = ABS( tr(ji,jj,jk,jn,Kbb) ) / ( ABS( tr(ji,jj,jk,jn,Krhs) ) + rtrn )
                  xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
               ENDIF
            END_3D
         END DO
         !                                ! where at least 1 tracer concentration becomes negative
         !                                ! 
         ! Concentrations are updated
         DO jn = jp_pcs0, jp_pcs1
           tr(:,:,:,jn,Kbb) = tr(:,:,:,jn,Kbb) + xnegtr(:,:,:) * tr(:,:,:,jn,Krhs)
         END DO

With the changes in #45 it would be fairly straight forward for me to add a (optional) similar check to the tracers integration where it multiplies all of the tendencies by xnegtr if we thought that was a good idea?

I'll have to check that this still conserves mass when individuals exert a tendency on a tracer as I imagine it may not (I think it would cause the system to gain mass unless the integration of individuals properties was added into the same sort of system, which it could and maybe should be), but this sort of thing could be a good idea anyway.

[iuryt]

Hi @jagoosw

Not sure this is going to be of any help, but I have a code snippet for that here.

This is also done in the biogeochemical modeling using PSOM (from Amala's group).

I noticed that, of course, if the gradients are extremely high by the initial conditions or something else, the machine-level negative values propagate and blow up everything and even zeroing negative values won't help.

In my simulations I checked and it was conserving mass (with some very little oscillations). I think that would be interesting to somehow inform user about mass conservation problems, or add some extra output for that?

[jagoosw]

Hi,

Thank you very much for this. It's useful to know that it conserved mass doing it that way. I'll have a play with something like that in LOBSTER and PISCES.

It might be useful to add a callback every so often that checks the mass and shows a warning if it's changed by more than a certain amount?

[johnryantaylor]

Thanks @jagoosw and @iuryt. Having this capability is a great idea. Negative tracers can cause problems in some cases and it would be great to at least have the option to maintain positive values. Regarding the mass/tracer conservation, its a good idea to have a diagnostic for this (like we do for nitrate in the LOBSTER model). I think that a diagnostic would be preferable to a warning though since there are likely going to be use cases where the tracers aren't conserved (e.g. adding a forcing term to mimic nutrient upwelling).

[glwagner]

Another approach is to use the "Patankar trick" described here:

https://www.sciencedirect.com/science/article/pii/S0168927403001016

This is not a strict guarantee of positivity, however, because it also requires both positive-preserving advection schemes and for diffusivities to play nice (if diffusion is treated implicitly).

[johnryantaylor]

Thanks @glwagner, I hadn't read this paper, but this looks like a very nice (and seemingly easy to implement) trick - definitely worth a try!

[jagoosw]

Thanks for this, it looks really interesting!

Since it needs the production and destruction terms separately, to avoid having to change lots of things like calculate_tendencys! perhaps instead we could implement it by having forcing functions return imaginary numbers and a new Timestepper that uses the components as P and D?

[glwagner]

I think we should add two terms to the RHS of the Oceananigans tracer equation: production and destruction, in addition to forcing. Implementations of BGC provided by this package can then use those slots to implement BGC terms, leaving the forcing slot open to users and their specific science applications. (It comes to mind that we may even want a third BGC-specific term for things like settling, when we can't identify whether it's production or dissipation.)

For specifying production and destruction, we may need to add a biogeochemistry property to Oceananigans models too.

We can then design a time-stepper + user interface for treating the destruction terms implicitly in time. This would be similar to the interface we currently have for the turbulence closure term (ie whether or not to use VerticallyImplicitTimeDiscretization).

This is probably not a small change to Oceananigans, since we might have to refactor / repurpose the vertically-implicit diffusion solver too, which is currently in Oceananigans.TurbulenceClosures:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl

I think though some kind of comprehensive design like this is probably necessary for long-term sustainability...

[johnryantaylor]

Hmm...I think that I spoke too soon when I said that it would be easy to implement :) I was thinking that we could have the user specify the production and destruction terms and then combine these into a forcing term, but now I realize that this approach won't generalize for different time-stepping methods.

[glwagner]

It's certainly straightforward to implement if you can change the time-stepper (which we can!)

But I think maybe you're right to say it's not "easy" because to do this well, we should design a generalizable interface in Oceananigans that OceanBioME can plug into.

Hard to say whether implicit time-stepping is preferred over the NEMO method. The NEMO method is obviously wrong, since the tendencies are distorted when tracers are close to zero. On the other hand, using an implicit time-stepping method doesn't "solve" the core issue (which is that the time-step is too long to resolve the process of tracer depletion when it's close to zero). The implicit method may be less computationally expensive though --- especially if we are already using vertically-implicit time discretization for diffusion.

I suspect that this kind of thing is only the beginning of a host of specialized numerical methods that might be useful for reacting tracers in Oceananigans.

[jagoosw]

In the mean time before we can implement a proper positivity preserving time stepper I did some experiments with different ways to force positivity. My test case is a system that I run for a year and where for around 100 days the ammonia is very low. I can make this run with no negatives if I set the time step to around 1.5minutes and am testing the methods with a 5minute time step (so this is probably a fairly extreme scenario and with smaller timesteps the errors would be smaller). My results are:

• Changing any negative value back to zero with a callback - system gains mass (as expected), over the year gained about 1.2% in total N and primary production was about 1.5% higher. The pattern of change to the system was also as you'd expect if there was slightly more ammonia available - slightly more growth and then slightly less death and particle production.
• NEMOs method of re-scaling the values of the other (relevant - in this case P, Z, NO $_3$, NH $_4$, D, DD, DOM) tracers in a cell to get negative tracer to zero and conserve budget. I implemented this by, in each cell, setting C $_f$ = 0 if C $_i<0$ otherwise C $_f$ = C $_i \cdot T/P$ where T is the sum of all tracers in the cell and P is the sum of positive tracers. This distributes the negative between all of the other tracers proportional to their original concentration, and since $\Sigma C_f = T P/P = T$ it conserves the total mass. This resulted in the same conservation of mass as the short time step, only increased total primary production by about 0.8%, and produced a much shorter lived, but occasionally more intense, system response.

This is the total budget as a fraction of the initial budget, the two small ones are $\sim10^{-11}$ random fluctuations but you can see the increase in budget of just resetting to zero.

Here is the primary production difference, you can see that when we force to zero its just like adding more nutrients, where as re scaling is more of a perturbation that it has an oscillatory (although still biased positive - perhaps todo with other factors like light availability, maybe experimenting with constant light would be better) response to:

And finally some poorly formatted heatmaps of the systems - original showing that the change when we rescale is much more contained but that there is a larger P and NO $_3$ response:

From this I tend to think that re-scaling is the better temporary solution @johnryantaylor @syou83syou83? We could interfere even more and only shuffle negative amounts between the sub-groups, e.g. if NH $_4$ is negative slightly decrease NO $_3$, if D is negative decrease DD, etc.?

[jagoosw]

To add to this, I've finally checked that advective forcing that smoothly goes to zero at the bottom: a) conserves the budget to the same degree as without (with a short time step)

b) By scaling negatives the budget is conserved as above and there is a similar response in primary production with no drastic changes in the whole system

Short timestep, no negative protection:

Longer timestep, scaled negatives:

Difference (scaled - no protection):

Sorry for all the messages, mainly adding this for the record!

Note for my own record: have now confirmed that these results are also reproduced for irregularly spaced grids including the velocity going to zero smoothly enough for stability.

[jagoosw]

Also, I think the reason this hasn't really been an issue before is because the time that NH4 gets close to zero is the same as when the mld is shallow so we've been dropping the time step for the diffusive CFL anyway

[glwagner]

Nice work! Interesting to see the effect of the different methods on dynamics. Can you also generate a "truth" run that uses a very very small time-step (small enough to resolve dynamics / prevent negativity)?

[jagoosw]

Thanks! This one is like that with a short time step that stays positive by its self. The differences are compared to runs with short time steps like this.

Short timestep, no negative protection:

I can plot the short time step one without advective forcing but the differences are almost identical (except at the bottom of the domain).

[johnryantaylor]

This is great, thanks @jagoosw! It might be useful to plot the difference as a percentage of the total, but they look quite small (a few percent which is arguably well within the 'tolerance' given many uncertainties in the BGC system).

[jagoosw]

Here are plots with no sinking as percentages:
Zeroing:

Scaling:

Because most of the difference in the NH $_4$ is when the concentration is very small the percentage in that particular cell can be very large. Here are the results as a percentage of the maximum reached over the whole time in each particular cell:
Zeroing:

Scaling: