Negative passive tracer concentration

[ilicakme]

Hi everybody,

I got a problem about passive tracer concentration.
I release some dye in my simulation as
t(i,j,k,1,itrc) = 100.0_r8

Dye is located as a box volume. It is 100 in a specific location
and zero everywhere. I was using 3rd order upwind scheme and I got negative values, so we thought maybe it is because of Gibbs oscillations so I used MPDATA and I didn't get any negative values.
After that I used point source rather than initial condition;
so initially dye is zero everywhere and I put some source terms as
Tsrc(is,k,itrc)=100.0_r8

I got two questions;
1) dye values in the simulations are very small order of 10^(-5). Why?
2) There are negative values in both MPDATA and 3rd order upwind scheme.

We thought maybe, when the dye is put as source, it is considered as flux, so if it is true, how can I plot only concentration?

Thanks,
Mehmet

[inga]

Hi,

I was about to post a problem about the negative concetrations as well. Few days ago I started to experiment with biology and dye and I discovered that putting as initial conditions 10.0 concentration of dye at some depth and 0.0 elsewhere I get negative concentrations down to -1.0. To discover the problem I looked into the roms code and at different stages (prestep,mixing, advection) put the max(0,tracer) conditions at the beginning and end of the procedures. I found out that it is mix4 that produces the negative values. It happens with all 3: t3dmix4_iso.h, t3dmix4_s.h and t3dmix4_geo.h procedures . With max(0,tracer) at the end of those schemes, no negative values are produced anymore.

Well, it seemed to me terrible that the horizontal mixing scheme produces negative values as it absolutely shouldn't, by all fluid dynamics principles.....

I tired to look deeply into the code but owing to my limited knowledge on programming I haven't yet succeed.

[wilkin]

I found out that it is mix4 that produces the negative values. It happens with all 3: t3dmix4_iso.h, t3dmix4_s.h and t3dmix4_geo.h procedures .

Horizontal mixing indeed should not violate monotonicity of conservative model tracer concentration. Those of us who are experienced ROMS users have not noticed this is our applications. In particular, it is essentially impossible for the horizontal mixing along s-surfaces to cause this. Something else is wrong.

In order for users of this forum to help guide you toward resolving the origin of the problem we need to have much more specific information on how you meant to pose the physical problem and how you configured it numerically (not just the model options, but a description of how you have implemented your physical initial/boundary conditions in the discrete model sense).

Here are several considerations that might be contributing to the model behaviours described in this thread:

1. Sharp fronts: Initial conditions with sharp spatial gradients are a severe test of the model algorithms, particularly advection. The Gibbs phenomenon (especially in 2nd-order schemes) noted in the posting by ilicakme is a well-known example. This effect is exacerbated by sharp fronts. Keep in mind that you can't expect the model to resolve processes at the grid scale, and certainly not a 1-grid-point wide front. So setting initial conditions that are abruptly discontinuous on the grid scale are almost doomed to failure. In our successful dye-tracing simulations we initialize a smooth Gaussian shaped patch with length scale of at least a few grid spacings - the evolution of the patch is modeled well.

2. MPDATA should reduce the loss of monotonicity that can still occur with 3rd-order upwind. However, at point sources all the advection schemes get messy and I think several of us would acknowledge that the properties of the respective advection schemes may not be fully respected at a point source location; our results for river sources entering coarsely resolved 'estuaries' indicate there is still progress to be made on how this is implemented. Interior to the model, away from point sources, MPDATA should perform well. MPDATA is one of the most recent enhancements to the code and has lately undergone some debugging, so if you have not recently updated I encourage you to do so before working further with MPDATA.

3. Global conservation:

With max(0,tracer) at the end of those schemes, no negative values are produced anymore.

This quick fix will eliminate unphysical concentrations less than zero but will cause a loss of global conservation (integral over the entire domain) of total tracer. The advection schemes are all implemented to maintain the total inventory of a conservative tracer, so that an initial patch of X units of conservative dye will be maintained even if during part of the simulation some concentrations dip below zero or peak above max initial value. This returns to my point 1 above: focus can't be on isolated values at the grid scale but should be on the resolved scales (at least a few grid scales). The biology codes go to great lengths to respect global conservation. In fasham.h a copy is made of the bio tracers ...

Code: Select all

line numbers:
368 	!  Extract biological variables from tracer arrays, place them into
369 	!  scratch arrays, and restrict their values to be positive definite.
370 	!  All the tracer at input have transport units: m Tunits.

376 	              Bio(i,k,indx)=MAX(t(i,j,k,nnew,indx),0.0_r8)*             &
377 	     &                      Hz_inv(i,k)

... concentration > 0 is enforced, the bio processes act on these fields, then the total change in concentration due to bio processes is added back to the original tracer values ...

Code: Select all

line numbers:
1171 	!-----------------------------------------------------------------------
1172 	!  Update global tracer variables

1187 	              t(i,j,k,nnew,indx)=MIN(t(i,j,k,nnew,indx),0.0_r8)+        &
1188 	     &                           Hz(i,j,k)*Bio(i,k,indx)

... which means some values might still be less than zero, but the global inventory is conserved. This way an analysis of the total budget of, say, nitrogen, can be computed over a large spatial domain and exact balances between initial, final, uptake, sinks and sources, can be achieved. If the actual tracer array values were arbitrarily bounded with max(0,tracer) statements these would effectively become a source term and physical balances would be lost.

[mizuta]

Just for your information, equations with a biharmonic term (which is used for mix4... option) e.g.
p_t = -p_xxxx with p=delta(x) at t=0
do have negative solutions. Such features of the biharmonic term are quite different from those of the harmonic term (e.g. Pedlosky & Chapman 1993, JPO and related papers).
Genta

[ilicakme]

Thanks for responding,

But my problem was more point source tracer concentration rather than initial conditions.

Negative concentration at initial condition is probably because of Gibbs oscillations. Since when I used MPDATA there were no negative values.

But when I used point sources, even using MPDATA didn't help, and I got negative values. Can anyone help me about this?

Another thing is I noticed in Roms 3.0, in all the input files there are new 4 turbulence closures mentioned; k-kl, k-epsilon, k-omega, and gen (generic). However in Roms 2.2, there are different closures; MY2.5, k-epsilon, k-omega, and k-tao.

In my research I'm comparing different closures, so is there a particular reason to remove k-tao from the model?

[jcwarner]

For the negative conc issue: There is a possibility for negative tracer concentrations (even with MPDATA) due to the formulation of the point sources. The point sources are applied as fluxes during the time stepping. If the user specifies a flux that exceeds the concentration values in a cell, then more of that material will be transported out than exists. This will cause negative values. Even with mpdata this can happen.

It is possible to prevent mpdata from going negative for the point sources. But this would entail allowing the point sources to be coded a different way. I know other developers have discussed these issues. I do not know how to set priorities for these things, but it is possible.

Beside the possibility of point source issues, i do not get any negative values with mpdata in any of my tests. Recently in v3.0 I made a change with MPDATA to use a more correct grid cell thickness in the corrector time step. Let me know if this causes any problems for anyone.

For GLS you can refer to :
Warner, J.C., Sherwood, C., Arango, H., and Signell, R. (2005). “Performance of Four Turbulence Closure Models Implemented Using a Generic Length Scale Method.” Ocean Modelling, v. 8/1-2, p. 81-113.

I did not test GLS with the k-tau settings. I had included them before because that method exists. But knowing that i did not include any tests for k-tau in our OM paper, I did not feel comfortable providing those values as possible coefficients in the ocean*.in files.

You, of course, can try any set that you like, but be mindful that these coefficients must be chosen carefully. You can read the plethora of Umlauf and Burchard papers to see what those guys have been up to recently.
ROMS 3.0 has GLS, MY25, KPP, BVF, and ana_vmix as turbulence closure options.

-john warner, jcwarner@usgs.gov