Can we have only passive tracers in rivers/point sources?

[lanerolle]

I am attempting to set-up an application where some point sources will be rivers (with transport, T, S) but others will be passive tracer sources only (with transport, passive tracers but no T and S). Now when I look at the river_flag option in ROMS :

river_flag =0 is all tracers off
river_flag =1 only temperature is on
river flag =2 only salinity is on
river flag =3 both temperature and salinity are on

Is there a river_flag value where only passive tracer is on? Or, should we modify the ROMS on our own (e.g. via Frc in Moduldes/mod_sources.F, etc.) to include this option?

[kate]

Have you updated lately? From an input file:

Code: Select all

!  LtracerSrc  Logical switches (T/F) to specify which tracer variables
!              to consider when the option TS_PSOURCE is activated. Only
!              NAT active tracers (temperature, salinity) and NPT inert
!              tracers need to be specified here:
!
!                LtracerSrc(itemp,ng)     for temperature (itemp=1)
!                LtracerSrc(isalt,ng)     for salinity    (isalt=2)
!                LtracerSrc(NAT+1,ng)     for inert tracer 1
!                ...                      ...
!                LtracerSrc(NAT+NPT,ng)   for inert tracer NPT
!
....
!              This logical switch REPLACES and ELIMINATES the need to have
!              or read the variable "river_flag(river)" in the input rivers
!              forcing NetCDF file:

[wilkin]

As Kate notes, the functionality you seek is already available. The LTracerSrc logic was introduced in December 2009. See the notes associated with the trac ticket:
https://www.myroms.org/projects/src/ticket/376

The change is backward compatible with old river forcing files, but introduces convenient new functionality. Notably, the handling of which rivers have which tracer sources (inert, biological, sediment etc) can be removed from the netcdf file itself: there is no need any more to meddle with the

Code: Select all

river_flag = ... 

values in the nc file when experimenting with adding different sources.

However, you do need to add the passive tracers to ALL rivers in the netcdf file even if you don't intend to activate them in ALL rivers. The logic here becomes obvious when you consider the netcdf file dimensions.

For example:

Code: Select all

dimensions:
	river = 7 ;
	s_rho = 36 ;
	river_time = UNLIMITED
variables:
double river_salt(river_time, s_rho, river) ;
	river_salt:long_name = "river runoff salinity" ;
	river_salt:units = "PSU" ;
	river_salt:field = "river_salt, scalar, series" ;
	river_salt:time = "river_time" ;
double river_dye_01(river_time, s_rho, river) ;
	river_dye_01:long_name = "river inert tracer" ;
	river_dye_01:units = "non-dimensional" ;
	river_dye_01:field = "river_dye_01, scalar, series" ;
	river_dye_01:time = "river_time" ;

All sources have the dimension "river" - there is really no other flexible way to easily handle this. But if you don't want dye_01 in rivers 3,4,5 for example, then you just turn off those logical flags and values in river_dye_01(:,:,[3 4 5]) are never used.

To the best of my knowledge you must provide values for the active tracers, temp and salt, for all rivers. You can't have a mass flux entering and not a corresponding temperature and salinity. But I guess you could try turning off the temperature sources and see what happens.

For salt this is easy because the river salt is zero. If you are adding Q m3/s of volume transport then it makes sense that the zero salinity entering the source cell dilutes the salinity there. Similarly for temperature - so you need some notion of the river temperature. If you don't want the river to change the water temperature, then you need to conjure up some time series that will have minimal impact. But you are implicitly saying that you know the river temperature to be the same as the water body it flows in to - only you don't know that. So this will be tricky to implement.

For other tracers (e.g. biology state variables, for which there are now corresponding LTracerSrc flags in the biology .in files) if you have no information about the concentrations then the rational thing to do is make them zero. The mass flux of zero concentration will dilute that state variable in the source cell, as it should.

[lanerolle]

I am extremely concerned by your reply. So to summarize, from what you wrote, if we have a 3D baroclinic simulation with :

1. N1 rivers treated as point sources
2. N2 point sources in the middle of the ocean with Np passive tracers

We have to ensure that ALL of the point sources (both rivers and those in the middle of the ocean - i.e. N1+N2 sources) need to have T, S and Np passive tracers. Now for the rivers, we can put the passive tracers to the background/zero values so that they do nothing to the computation. However, for the N2 point sources in the middle of the ocean, in addition to the Np passive tracer values, we ALSO have to specify T and S for them!!!

How do we know the T and S values for these locations in the middle of the ocean in advance? If on the otherhand we do specify some reference T, S values for them, these values will come into the computation and chance the T and S fields in the vicinity of the point sources which in turn will change the vertical eddy-viscosity/stratification and hence modify the dispersion of the passive tracers (and also for T and S themselves)!!!

So in short, how do we configure the new ROMS to simulate a dye experiment in the middle of the ocean where rather than doing a one-time dump of dye, we have put in a dye source which emits dye at a say constant rate of Q m3/s? Alternatively, how can we use the ROMS to do say an deep sea oil spill where oil (to be treated as a passive tracer) is being spilt in to the ocean (at a particular depth say) at a constant rate of Qx m3/s? Another example would be the calculation of say Phosphate residence times/budgets where some Phosphate will enter the ocean via rivers and the other via point sources (where only Phosphate time-series observations/measurements are avilable and not any T and S time-series)? Will we also need to know some T, S values for these locations in advance to do these simulations?

Please let me know what you think.

[wilkin]

Like I said:

I guess you could try turning off the temperature sources and see what happens.

Run a simple test for an internal source like a well head leak with only one dye and see what happens. Is temperature conserved globally, or does the volume flux effectively dilute the heat as if the source had temperature zero? You could do this with a slightly enlarged BIOTOY case.

Let us know the results.

[kate]

I assumed he was saying that for sources for which you've turned off the T,S sourceness, you still need to provide bogus values of T,S in the file.

As for oil, a buoyant plume isn't going to be exactly passive. I saw a diagram with various amounts of mixing going on, resulting in sheets of oil blends at different depths as each mixture rose to its own density. Messy.

[lanerolle]

OK, I'll post the results of the numerical experiment(s) you suggested (in the next few days).

Meanwhile, I was thinking, can we do the following two changes to ROMS?

(1) change ocean.in so that LTracerSrc can take the following form? Assume that we have 5 point sources and we have 2 passive tracers. For the first river we have T, S only (no passive tracers), for the second river we have T only, for the third river S only and for the fourth river the first passive tracer only and for the fifth river the two passive tracers only. Then, we could change LTracerSrc in ocean.in to take the following form where the rows are the number of point sources (=5) and the columns are logical switches for temperature, salinity, passive tracer 1 and passive tracer 2.

The number of rows = number of point sources and the number of colums = 2 + number of passive tracers, where the 2 is the T and S - i.e. the two active tracers. So we have in ocean.in :


LTracerSrc = T T F F
T F F F
F T F F
F F T F
F F T T

What do you think? (I tried to indent the last four rows for clarity but it does not seem to work in the myroms.org editor)

(2) Alternatively, we can change the ROMS coding associated with the TS_PSOURCE CPP option where rather than looping over the range 1:(NAT+NPT) [i.e. the total number of active and passive tracers], we loop over 1:NAT and 1:NPT separately and only do one or the other or both depending on some switch within ocean.in or a new CPP option. I think lumping the active and passive tracers together and looping over them seamlessly is too restrictive.

[wilkin]

I don't think you're fully considering the logical consequences of a point source that introduces some Q m3/s volume flow into a grid cell.

For the first river we have T, S only (no passive tracers), for the second river we have T only, for the third river S only

Say your "first river" adds a flow of Q1 into the source cell. (This is represented as a volume flux through the face of the cell, and flow direction, that you designate defines the source). The volume of the water column is increased by Q*dt distributed vertically according to Vshape (less what flows out). If there is no dye_01 (say) associated with river 1, then I believe the volume flux dilutes the concentration of dye_01 (and all other tracers) that already exists in the cell.

You seem to be asking that a river source be allowed not to alter the concentration of a dye. The river source notion here in ROMS is not the same as adding a delta function to the RHS of the passive tracer equation. That would add tracer without affecting the model volume.

[lanerolle]

I did some ROMS simulations with a single point source. In the river/point source forcing NetCDF file, in the vertical, I put the distribution to be uniform and as I was only interested in a passive tracer evolution, I put both T and S to zero (in time and space). As I was also trying to simulate a "dump", I put the river transport to zero (in reality it should be a hat-function with a width of ~1-3 hours). The only time-varying quantity was the passive tracer (river_dye_01) - see the plot below. My background passive tracer value was zero and so was my open ocean boundary value. I edited the toms.in file so that I had :

LtracerSrc == F F T ! temperature, salinity, inert

The CPP options I used included :

#define T_PASSIVE
#ifdef T_PASSIVE
# define ANA_BPFLUX
# define ANA_SPFLUX
#endif

and, I tried two configurations : (a) #undef TS_PSOURCE and (b) #define TS_PSOURCE.

For (a), I was expecting nothing to happen and for the passive tracer to be constantly held at the background/open boundary value of zero. For (b), I was expecting something like a diffusive process to take place. The results are given in the plots below for the first 16-days of the simulation - i.e. days 68-84 (please note the difference in color scales for the two Hovmoller plots).

As you can see, even with (a), the passive tracer begins to change and evolve with time on day 78 (which is when the river/point source has non-zero passive tracer values). With (b) the situation goes completely out of control and we get huge positive and negative passive tracer values. Although the passive tracer from the river/point source forcing file has only positive values, the plots below show that they in reality take both positive and negative values.

Could someone explain to me what is happening in the ROMS and why we get positive and negative passive tracer values and why a non-zero passive tracer evolution takes place for scenario (a)?

Passive tracer time-history from the rive forcng file

Hovmoller plot of passive tracer with #undef TS_PSOURCE

[lanerolle]

Sorry, I forgot to add that in the simulations pertaining to (a) and (b), the resulting ROMS T and S fields between them are identical. So differences between them are only seen in the passive tracer field.

[m.hadfield]

I recently needed to introduce passive tracer sources into the interior of a ROMS simulation. I couldn't get the point source system to do what I wanted in a toy simulation (it seems to generate a dipole when the point source is not adjacent to the land mask) so I added appropriate terms to the RHS of the tracer tendency equation in step3d_t.F by copying and adapting the nudging code. I neglected the very small volume flux associated with my source. (Ideally one could use Q_PSOURCE for that, but my experiments confirmed what I had read on the forum, namely that it is broken.)

Mind you, based on my recent performance on the issue of float-restarts, I may not be the most reliable source of information on what ROMS does and does not do. As always, you need to do your own tests. I recommend small-domain test cases with forcing specified analytically.

[wilkin]

Lyon,

You say you use a uniform distribution in the vertical, by which I assume you are refering to the dye concentration profile.

But you didn't say whether you have UV_PSOURCE defined, or not, in both cases (a) and (b). Assuming you did, what was the vertical shape of the velocity profile?

Aside: Have you used Q_SOURCE at all?

John.

[lanerolle]

To answer your questions :

1. Yes, the dye concentration is distributed uniformely in the vertical across all of the sigma layers (at each time snapshot) in the river/point source forcing NetCDF file - e.g. in the river_dye_01(Nt, Ns, Nz) variable, the dye value at each time is distributed evenly in the Nz dimension. Furthermore, I am also using a uniform river_Vshape(Nz) function and it has a form like f(z)=1/Nz so that Sum(f(z)) = 1.

2. As I was attempting to simulate an instantaneous dye "dump", I decided not to use UV_PSOURCE. Furthermore, in the river/point source forcing NetCDF file, I put river_temp(), river_salt() and river_transport() all to zero for all time. What I was trying to achieve was for the dye to naturally "diffusive" in ROMS becuase although its volume is negligible, it is highly potent and was trying to simulate something like putting a drop of ink in a glass of water which has been gently stirred already.

3. No I have not used Q_SOURCE in this simulation. Should I be using it?

[lanerolle]

So I repeated the above mentioned simulation with the only differences being :

(1) I included UV_PSOURCE in the CPP options when compiling the code and,
(2) for river_transport, I used 1 m3/s - a small value.

Please find attached the outcome from ROMS. Now, we find that the passive tracer is much closer to zero (the default/background value) but it still goes both positive and negative.

Could someone please let me know how I can include a point source/river in my ROMS simulations which only generates a passive tracer and not any T or S? This used to be so easy to do in the older ROMS like version 3.4, 3.6, etc.!

Also, do you think there is a bug in the way the rivers have been coded up in the newer ROMS and should I write something under "ROMS bugs"?

[m.hadfield]

Could someone please let me know how I can include a point source/river in my ROMS simulations which only generates a passive tracer and not any T or S? This used to be so easy to do in the older ROMS like version 3.4, 3.6, etc.!

To implement a point source for selected tracers, use LtracerSrc in the input file.

Also, do you think there is a bug in the way the rivers have been coded up in the newer ROMS and should I write something under "ROMS bugs"?

I don't know if it's a bug or a feature, but as far as I can see point tracer sources (TS_PSOURCE) are useful only at points adjacent to the land mask or the edge of the domain. Elsewhere they generate a dipole, as I pointed out above (8 september). And point mass sources (Q_PSOURCE) are broken.

[lanerolle]

Yes, I use the LtracerSrc in toms.in and have put it to F F T so that I only activate one (single) passive tracer.

Yes, in the past, I have always put point sources adjacent to land. For this application, it is in the middle of the ocean and just like you, I too get a dipole with positive-negative values. So in a way, I find it comforting that we both get the same outcomes. However, now the problem is how to fix this and allow for point sources in the middle of the ocean (e.g. to simulate a broken pipe underwater which is pumping contaminants in to the ocean, etc.).

Thanks for letting me know that Q_SOURCE is broken - I was going to use it to see what happens and now I know that I should not do so.

[m.hadfield]

However, now the problem is how to fix this and allow for point sources in the middle of the ocean (e.g. to simulate a broken pipe underwater which is pumping contaminants in to the ocean, etc.).

My post of 8 September (which did not seem to draw any reaction at the time) refers to this issue:

I recently needed to introduce passive tracer sources into the interior of a ROMS simulation. I couldn't get the point source system to do what I wanted in a toy simulation (it seems to generate a dipole when the point source is not adjacent to the land mask) so I added appropriate terms to the RHS of the tracer tendency equation in step3d_t.F by copying and adapting the nudging code. I neglected the very small volume flux associated with my source. (Ideally one could use Q_PSOURCE for that, but my experiments confirmed what I had read on the forum, namely that it is broken.)

Adding terms to the RHS of the tracer tendency equation is not that hard: you have to get the scaling right, obviously. I have done it in a very ad-hoc way, which would break if I changed the horizontal or vertical grid spacing. The source strength is constant: time-variation opens a whole new can of worms, but as I see it one could adapt the existing point-source infrastructure and extend it with a new source type and a new preprocessor symbol (M_PSOURCE?) Hernan suggested I add it to the wish list. I will do that right now.

As I see it, your choices are:

• Write some ad hoc code yourself.
• Wait for someone (e.g. me) to write the M_PSOURCE code. This may well happen in the not-too-distant future. Alternatively it may not.
• Write this more general code yourself.

[wilkin]

As Mark points out, the point source option is coded for the situation where the source lies on a land/sea boundary. Specifically, it causes the mass (or volume) flow and tracer to enter through a velocity face (u, or v) of a rho-centered cell.

If you look at the code, you'll see it achieves this by modifying the normal velocity on the cell face (which would otherwise have been zero because of the land/sea mask boundary condition) to take the value of the river velocity. This is in step3d_uv.F.

The influence of tracers at the source is introduced in step3d_t.F. If you look at the code you'll see ROMS takes the velocity at the cell face (which would have been modified in step3d_uv to be the river flow) and computes the tracer flux through that cell face. This is Huon*Tsrc (Huon is the vertical grid metric Hz, times velocity u, divided (o for over) n (the reciprocal horizontal grid metric "dy"). Just think of Huon as the velocity in flux form.

If you configure ROMS to place the source in the middle of the ocean, and not on the coast, and do not activate the velocity source, then Huon is simply the velocity that ROMS has computed for that time step. Multiplying by the tracer value (Tsrc) you will get fluxes on that face whose sign will depend on u. Obviously, this makes no physical sense and does nothing like what you seem to think this configuration would.

This is not a bug, or a feature. You are trying to accomplish something for which this code was never designed and cannot achieve.

If you want to add a right-hand-side forcing term - a "source" of magnitude Q - to the tracer equation that is effectively a delta function in space, then it needs to be added at the cell center (logically a "rho" point). You'll have to introduce this in step3d_t.F somewhere near where the TCLM_NUDGING is done. Something like t(i,j, ... nnew) = t(i,j,...nnew) + Q.dt but scaled by the appropriate metrics. That will depend on how you define Q.

But even this seems ill conceived - what sort of physical source to you envision acts this way in the middle of the water column?

A more likely physical scenario is that there is a source at the seafloor or sea surface - like a local geothermal source of heat. In this case, you would want to impose a bottom flux of some magnitude. ROMS provides for this by allowing you to see bottom passive tracer fluxes. By default these are set to zero everywhere in analytical.F, but you can customize the ana_btflux.h to do so.

[m.hadfield]

But even this seems ill conceived - what sort of physical source to you envision acts this way in the middle of the water column?

Usually sources like this are intended to represent the footprint on the model-resolved scale of some subgrid process. A wastewater pipe? A cloud of sediment stirred up by dredging? A fish farm releasing nutrients into the upper water column? Sometimes there are significant unresolved processes that need to be resolved before the model-scale source is constructed (e.g. buoyant plume effects), sometimes not so much.

[wilkin]

Mark,

All of those scenarios would best be represented - in ROMS - as a rho-point-centered source of tracer. Not as a horizontal flux at the vertical face of a grid cell in the middle of the water column.

Several of those scenarios could very sensibly be represented as surface or bottom fluxes of tracer to the vertical diffusion term. This option already exists.

Evidently, we need to fix the Q_PSOURCE option.

John.

[m.hadfield]

Mark,

All of those scenarios would best be represented - in ROMS - as a rho-point-centered source of tracer. Not as a horizontal flux at the vertical face of a grid cell in the middle of the water column.

Several of those scenarios could very sensibly be represented as surface or bottom fluxes of tracer to the vertical diffusion term. This option already exists.

Evidently, we need to fix the Q_PSOURCE option.

John.

A rho-point-centred source of tracer is exactly what I am proposing.

I thought Q_PSOURCE was intended to provide a volume source only, not a tracer source. A quick look through Nonlinear/step3d_t.F has cleared up my confusion. Yes, we need to fix Q_PSOURCE.

As you say, a surface or bottom flux may be an option also: it depends on the specific situation.

[arango]

So what is wrong with Q_PSOURCE? This is not in my list of ROMS problems. This type of source is complicated because involves the continuity equations of ROMS. I am always uneasy messing around with the 2D and 3D continuity equations. It may explains why this is broken. This is very delicate stuff.

[m.hadfield]

So what is wrong with Q_PSOURCE? This is not in my list of ROMS problems. This type of source is complicated because involves the continuity equations of ROMS. I am always uneasy messing around with the 2D and 3D continuity equations. It may explains why this is broken. This is very delicate stuff.

I did a test case, a small enclosed domain, with a single volume source (Q_PSOURCE) at the centre. The volume in the domain increases at the required rate, but there is a column of negative vertical velocities at the source location, which seems wrong. I haven't checked the combination of TS_PSOURCE and Q_PSOURCE. I will have another look.

[johnluick]

Mark, I had the same problem until I restricted the inflow to the top bin or few bins (I forget just now). I would be interested to know the outcome if you were to try this (top bin only).

But what if inflow is at or near the bottom, like the brine from a desal plant. Has someone sorted this out?

[arango]

Every tracer in ROMS needs surface and bottom boundary conditions: stflx(x,y,itrc) and btflx(x,y,itrc). The bottom tracer flux is usually set to zero. We use the ANA_BTFLUX CPP option to assign such values. See routine ROMS/Functionals/ana_btflux.h. These fluxes are used as boundary conditions to the vertical diffusion equation. The surface flux for temperature is the surface net heat flux. For salinity we have surface freshwater flux. Other passive tracer have surface and bottom fluxes with units m/s times tracer concentration.

Similarly, vertical boundary conditions are needed for the 2D and 3D momentum equations in terms of surface and bottom stress: sustr(x,y), svstr(x,y), bustr(x,y), and bvstr(x,y).

Fluxes at the surface and bottom are the more straight way to force any fluid. Sometime we can use these fluxes as a simple body force to avoid all the issues involved with a vertical mixing parametrization.

[m.hadfield]

Mark, I had the same problem until I restricted the inflow to the top bin or few bins (I forget just now). I would be interested to know the outcome if you were to try this (top bin only).

But what if inflow is at or near the bottom, like the brine from a desal plant. Has someone sorted this out?

I've realised that I may be confused about what Q_PSOURCE is trying to do. I thought a positive value implies an upwards volume flux (and so a positive Q producing a negative vertical velocity seems wrong), but perhaps it's supposed to be a downwards flux??? And in any case, you won't get a tracer source term (with Q_PSOURCE and TS_PSOURCE) unless you have a divergence in the volume flux. So I need to do some tests. This will take a day or 2 as I am currently at home nursing a cold.

[lanerolle]

So taking Mark's and John's advice, I set-up a passive tracer simulation where the point source was adjacent to land (to avoid the horizontal dipole). Now, I have a new problem!!! The passive tracer assumes both positive and negative values in the horizontal, vertical and also in time (intermittently)! I used the following CPP options (in addition to the usual):

#define T_PASSIVE
#ifdef T_PASSIVE
# define ANA_BPFLUX
# define ANA_SPFLUX
#endif

#define NONLIN_EOS

#define UV_ADV
#define UV_SADVECTION

#define TS_SVADVECTION
#define TS_U3HADVECTION

#define DJ_GRADPS
#define SPLINES

#define ANA_BTFLUX

#define TS_PSOURCE
#define UV_PSOURCE

I thought because we use an upstream-biased advection, we do not get oscillatory tracers and momentum in space. Also because of the Adams-Moulton-Leapfrog predictor corrector time stepping, the parasitic/computational/numerical temporal mode is damped and hence we should not have +/- values in time either. But this does not appear to be the case.

I set-up a series of simulations using the passive tracer time-series shown in the Figure 1 and they are with :

(a) river transport = 1000 m3/s, Vshape = 1/Nz, P_amb=0.0
(b) river transport = 1000 m3/s, Vshape = [zw(k)-zw(k-1)]/H (sums to 1.0), P_amb=0.0
(c) as in (b) but with river transport = 50 m3/s, P_amb=0.0
(d) as in (c) but with P_amb=1.0e-6

where P_amb is the value employed as the initial and open boundary passive tracer level and the minimum value for the point source time series also. The river transport was ramped-up linearly in time over a period of 1-day.

The results are shown in Figures 2-5. In Figures 2-4, I have plotted a Hovmoller plot (depth-time) at two grid points East of the point source (the river transport > 0 and river direction = 0 so the point source flow is to the East). Figure 5 is the result from run (c) and is a horizontal spatial plot of the vertically integrated passive tracer at each grid point. In these plots, I have artificially changed the color scales to highlight the positive and negative values. Using P_amb=1.0e-6 did not help to alleviate the negative passive tracer values and run (c) and run (d) gave practically the same outcomes.

So as you can see, ROMS generates oscillatory positive/negative tracer values in space and time and using different background tracer values does not help. Using different Vshape functions does not help either. Furthermore, it appears that for smaller river transports, we get more negative values than for larger transports. For my application, ideally, I need to use a very small river transport value (prescribed by observations and it is much less than 1000 m3/s).

Could someone please let me know why this is happening in ROMS and suggest some workarounds for me? I need to do this simulation urgently and it is high priority.

[arango]

I have told you several times in the past to not use TS_SVADVECTION and UV_SVADVECTION in realistic applications I mentioned this when we visit NOAA and every time that the subject comes up. I see that you have ignored the advice. This vertical advection options are in ROMS legacy code and should be only used in highly idealized applications Perhaps, I should remove these options from the code to avoid repeating myself all the time.

I have also recommended in all ROMS workshops to use TS_A4HADVECTION and TS_A4VADVECTION with explicit diffusivity (TS_DIF2 and MIX_TS_GEO). For vertical advection of momentum you should use UV_C4VADVECTION.

Just try the recommendation and compare the solutions. Then, you will know what I am talking about. This is the best advice that I can give to any ROMS user...

[lanerolle]

OK, I'll try these CPP options and post the results I fine. It appears that the CPP option you mentioned - UV_C4VADVECTION is not used in ROMS! It is documented in Utility/checkdefs.h but it does not appear anywhere else in the code and its not written-up in Include/cppdefs.h either. I am using ROMS svn version 562.

[arango]

Yes, it is the default option. I like to add it to remind me what I am using for vertical advection of momentum.

[lanerolle]

I have a quick question : the "optimal" CPP options you have recommended are 4th order accurate in space. Therefore, to remove the spatial oscillations, shouldn't we use a biharmonic viscosity/diffusivity (e.g. TS_DIFF4, etc.) instead of the harmonic options (e.g. TS_DIFF2) you suggested? - i.e. are 2nd order harmonic viscosity/diffusivity operators able to filter out oscillations of 4th order numerical methods?

[lanerolle]

So I set-up a simulation using exactly the CPP options which Hernan recommended. My simulation is such that dt=15s (baroclinic), dx=dy=1000m (approx). For the diffusivity, I tried K=3, 7, 67, 100 m2/s which correspond to the dimensionless parameter mu=K*dt/(dx*dx)=4.5e-5, 1.05e-4, 1.0e-3, 1.5e-3 all of which are quite small specially relative to the CFL number which is typically ~0.5.

In all cases, the model blows-up almost immediately! I tried using K > 0 only for T with K=0 for S and passive tracer (P), using K > 0 for S only with K=0 for T and P and having K > 0 for P only with K=0 for T and S. In all cases the fields quickly went NaN. I then tried cutting dt by 1/3 so that dt=5s and ROMS ran a little longer but still blew-up.

So it appears that the simulation is unconditionally unstable! For this domain, the bathymetry spans 5m to 5500m (very deep and very steep) and there are also many naturally occuring tall, steep seamounts.

Should I be using MIX_S_TS instead of MIX_GEO_TS due to the steep topography and seamounts? Also, I saw that there is a CPP option called MIX_STABILITY in t3dmix_geo.h but it has been undefined for some reason - is it worth using it?

Please let me know what you think. Any and all advice is most welcome.

[lanerolle]

Sorry, I should have mentioned that with K=0 for T, S and P, ROMS runs successfully in a stable fashion. However the P field has a lot of + and - values - even more so that using the upstream-biased advection CPP options. ROMS only blows up when I used K > 0 - even small values such as 3 m2/s.

[arango]

It is possible that you have a vertical CFL violation due to the rotated diffusion tensor. The constraint on time-step is higher with the rotated tensors when dealing with steep and tall bathymetry. It even get worse with the biharmonic operators. Implicit time-stepping of these algorithms, in both the horizontal and vertical discrete grids, is very tricky in parallel applications. I want to look at this issue in the future...

We need to use the rotated diffusion tensor in steep and tall bathymetry. Otherwise, you will have excessive numerical diapycnal diffusion that tends to over mix water masses. We need to take a look at the evolution of T-S properties in deep waters. This problem is inherent to all terrain-following coordinates models. I have investigated this in the North Atlantic. Others, have looked at this problem in the North and South Pacific. There is literature about this subject and one has to be aware it when configuring this type of application. I have showed this at several ROMS Workshops, every time that I have new results.

In my experiments, I get the best behavior when I use TS_A4HADVECTION, TS_A4VADVECTION, TS_DIF2, UV_VIS2, MIX_TS_GEO, and MIX_S_UV. I also use Vtransform=2, Vstretching=2 or 4, theta_s=7.0, theta_b=0.1 or 0.2, and Tcline between 250 to 500 m. Sometimes, I need to work smoothing the bathymetry a little with in the hot spots. The third-order upstream bias advection (TS_U3HADVECTION) for tracer may lead to a lot of implicit diapycnal diffusionin in some applications and we cannot control it. We just need to change the tracers advection scheme. The 4th-order Akima scheme doesn't have this but we need to use both viscosity and diffusion.

Now, the viscosity and diffusion in ROMS are implemented as physical processes of the fluid and not just a numerical noise reduction. The viscosity is implemented in terms of the deviatory stress tensor whereas the diffusion is implemented in terms of the Redi tensor (See Wajsowicz (1993), Sadourny and Maynard (1997), Griffies and Hallberg (2000), and others). The diffusion and stresses exerted on the ocean fluid are real. So you need to look this from that point of view instead of just a numerical stability or noise reduction (smoothing) issue.

Your application is somewhere in Japan and I don't know well the bathymetry of this area neither I have set-up regional applications there before. I don't have any expertise in this area to give you any advice. In a situation like this, I always check the literature before setting any application to learn the experience of others.

Setting a complex application takes time, requires curiosity, and careful testing of the solutions to insure realistic representation and accuracy. There is a lot of issues that we need to check systematically. The results that you are showing above indicates to me that this application is not properly configured yet. You need to start looking at your configuration systematically since there are many parameters and issues to consider. Every application is unique and we need to figure out the proper configuration parameters and flags. What it works in some areas may not work or be relevant in others.

Good luck...

[lanerolle]

I have done an extensive series of experiments and here is what I found out (at least for my application) :

1. If I use MIX_GEO_TS with TS_DIF2 (or TS_DIF4), ROMS is unconditionally unstable! I tried using diffusivity values 0 <= K <= 100 m2/s and with time steps dt, dt/3, dt/5 and in all cases, ROMS blows-up almost immediately (first 3-10 time steps either from a cold start or a hot start/restart) with the exception of K=0. For K=0, the model is stable for dt, dt/3, dt/5 and runs indefinitely but the passive tracer has +/- values both in time and space and are so bad that the numerical solutions are practically unacceptable. For my application, dt=15s and NDTFAST=60 (I also tried values 20-60 and get the same result). I also tried using the experimental ROMS CPP option MIX_STABILITY but it did not help either.

2. Out of sheer desperation, I used MIX_S_TS with TS_DIF2 and then repeated my runs for 0 <= K <= 100 m2/s with dt=15s. In all cases, the model is stable and does NOT blow up!!! I do not need to cut down dt to get the simulations to be stable. Furthermore, I have looked at the results and the diffusivity acts in the expected manner - the more diffusivity the more smooth the solutions are (even if they are unphysical) and the less smoothing the more osciallatory (in space and time) they are.

Please see the attached plot. The grey line is the coastline. I have plotted the passive tracer field (at the surface) using a Log10 scale for improved clarity and when tracer values are negative, I set them to NaN before taking the Log (in Matlab). So negative values (to the right of the grey line) appear in the plot as white pixels and non-negative (i.e. positive) values have various colors. As you can see, the negative values disappear with increasing K and I have confirmed this to be true by examing the full 3D x time passive tracer fields. [In the plot, I chose K=0, 6.667, 33.333, 66.667 m2/s because my dt=15s, dx=dy=1000m (approx) and so mu=K*dt/(dx*dx)=0, 1 x 10^-4, 5 x 10^-4, 1 x 10^-3]

So it appears that there is a problem with the Geopotential horizontal mixing in ROMS and it renders the whole kernel unstable (at least for my application were there is realistic T, S stratification). On the otherhand, the classical Sigma-surface mixing works well and in the expected manner but we know that it adds an artificially large amount of smoothing.

Could this possibly be a formulation problem or a bug in the ROMS Geopotential horizontal mixing algorithms?

[johnluick]

Interesting - but mine is stable with TS_DIF2, MIX_S_UV, and MIX_GEO_TS. Same as yours with the addition of MIX_S_UV (unless you have that and didn't mention it). I just mention it in case it helps you track down the problem. John

[simion1232006]

I used MIX_GEO_TS and have not found any problems yet. Maybe it is case by case.

[wilkin]

I routinely use MIX_GEO_TS - it's fine. That is not your problem.

How are you controlling nonlinear instability in momentum? With MIX_S_UV or a horizontal momentum advection scheme that dissipates the short wavelengths?

[lanerolle]

As I was mainly interested in simulating a passive tracer, I wanted to minimize the "damage" I would do to the numerical solutions by introducing "numerical/articificial" viscosity and diffusivity. So what I did was :

1. For tracers, I used TS_DIF2 and MIX_S_TS (as MIX_GEO_TS gave blow-ups) and used diffusivity only for the (single) passive tracer - hence, in ocean.in, I had TNU2 == 0.d0 0.0d 33.333d0 (TNU4 == 3*0.d0)

2. For momentum, I did nothing - i.e. I did NOT define UV_DIF2 or MIX_GEO_UV or MIX_S_UV or any such option.

I also based the above on the seamount test problem which uses both MIX_GEO_TS and MIX_S_UV but then sets the viscosity and diffusivity coefficients in ocean_seamout.in to 0.d0 in all cases. Why doesn't this problem use MIX_GEO_UV instead of MIX_S_UV? Is there a reason for that?

Furthermore, Hernan has always advised me that the way ROMS does the diffusivity for tracers and viscosity for momentum is such that the coefficients do have a physical meaning and we need to choose realistic values. I do not know what the real values for these coefficients are. As for the passive tracer, I merely want to use a diffusivity to ensure that it is positivity-preserving and hence I use diffusivity in a purely numerical manner.

So do you advice me to use the following in ROMS?

#define TS_DIF2
#define MIX_GEO_TS
#define UV_DIF2
#define MIX_GEO_UV

(my advection CPP options are exactly as Hernan mentioned in his earlier reply in this thread)

and if so, how should I choose appropriate values for TNU2 and VISC2?

[arango]

The horizontal diffusion tensor is rotated along geopotentials or isopycnals. This in turns puts a time-step limitation due to the density slope approximation of the rotated tensor. The numerical scheme is not implicit so you need to lower the time step accordingly to satisfy the CFL (usually vertical) condition. This not necessarily implies that the system is unconditionally stable It is just that you need to use a very small time-step due to you stratification and flow regime. Therefore, you need to be more precise when using such terminology.

I have been using rotated (geopotentials and isopycnic) tensors in all my realistic applications for years without problems with both harmonic and biharmonic. I just need to use the appropriate time-step and viscosity with the stress tensor.

By the way, when doing isopycnic rotation we first rotate the tensor along geopotentials so the time-step limitation is even more stringent.

[lanerolle]

I have been using rotated (geopotentials and isopycnic) tensors in all my realistic applications for years without problems with both harmonic and biharmonic. I just need to use the appropriate time-step and viscosity with the stress tensor.

So In addition to MIX_GEO_TS + TS_DIF2 for the tracers, should we also use some value of viscosity for the momentum equations (as also John may have suggested)? If so :

(1) What are the best CPP options - MIX_GEO_UV + UV_DIF2? or MIX_S_UV + UV_DIF2?
(2) Why is MIX_S_UV CPP option sometimes preferred (e.g. Seamount problem) over MIX_GEO_UV? and,
(3) How should we pick a value for VISC2 (in ocean.in)? Based on numerical considerations (e.g. smoothness of solutions, stability, etc.) or based purely on physical considerations? If the latter, it may be quite tricky because we will need peer reviewed papers/observations for each and every region of the ocean around the world and VISC2 will also vary depending on seasons and meteorological conditions.