passive tracer not conserved with TS_MPDATA in shallow water

[schen]

Hi all,
I set up a very simple case to test whether passive tracer is conserved in shallow waters. I used sediment as a tracer (only one size class) but set the setting velocity to zero. So it behaves like passive tracer. When I used the ROMS-default advection scheme, the total mass conserves very well (see figure below). But, when I used TS_MPDATA, the total mass is not conserved.

Setup:
In Fig 1(a) below. A beach with a constant slope (slope is 4m/6km). The north, south, and east boundaries are closed. At the western boundary, I imposed M2 tides with WEST_FSCLAMPED + WEST_M2REDUCED + WEST_M3GRADIENT (tidal amplitude is 1 meter). There are 243 x 48 x 10 grid points (east-west, north-south, along-shore, vertical, respectively). The dx and dy are 100 and 200 meter, respectively (x is west-east). Maximum water depth is 13 meter.

Initial condition:
In Fig. 1(b) below. There is an uniform patch of tracer (concentration is 0.1 kg/m3) from grid point i = 100 ~ 110, and j = 1~ 48. So, this is a 2D problem

Results:

Fig. 1(c), total mass conserves very well when the default advection is used.

Fig. 1(d), mass is NOT conserved with TS_MPDATA

Sensitivities:

I changed both barotropic and baroclinic time steps, but the results remain unchanged
DT = 60, NDTFAST = 20
DT = 60, NDTFAST = 60
DT = 30, NDTFAST = 30
DT = 10, NDTFAST = 10

Any suggestion is highly appreciated!!

[schen]

re-post the figure

[linzhenhua]

Hi,

I do not know the reason of the different results you get.But I have a question,how can you anticipate tracer to be constant with open boundary used?

And personally I do not think the MPDATA result is inferior to Default result,at least it could reflect the tidal signal,and on average the two results are nearly the same.

[schen]

I forgot to mention that the patch of suspended sediment never reaches the open boundary within the first two days.

John Warner has helped identify a potential bug in sed_settling.F
He suggested:

line 177 (or close to that number):

Code: Select all

              qc(i,k)=t(i,j,k,nnew,indx)*Hz_inv(i,k)

relace that line with:

Code: Select all

# ifdef TS_MPDATA
              qc(i,k)=t(i,j,k,3,indx)
# else
              qc(i,k)=t(i,j,k,nnew,indx)*Hz_inv(i,k)
# endif

I recompile, and it works!!

[arango]

Yes, good cash. Thank you. I updated the repository.

[pcheng]

I guess there is a bug in the diagnostics module in step3d_t.F. It forget to add horizontal advection term in diaTwrk when TS_MPDATA is defined.

After line 715, it may need to add:

# ifdef DIAGNOSTICS_TS
DiaTwrk(i,j,k,itrc,iThadv)=-cff1
# endif

Please take a check.