Negative suspended load concentrations near river runoff

[c.drinkorn]

Hi,

I get very low negative suspended load concentration values (down to several hundreds of kg/m^3) near river runoff. I have no sediment input via rivers. But the initial concentration in the water column is zero so I wonder what exactly a negative concentration should mean. Is there a way to prevent this? So far I just cut off anything below zero in the fields in order to make the "real" concentration visible. But I don't understand how these negative values come to be...

Thankful for any hint!

[jcwarner]

There are many reasons why there can be negative concentration values, but most likely it may be from the advection scheme. Suggest you use TS_MPDATA, that should remain positive. Also, we have added a new scheme called TS_HSIMT into COAWST that is also positive definite, it is based on some flux limter approaches (see Wu and Zhu https://doi.org/10.1016/j.ocemod.2009.12.001).

Most importantly though: you need to add up all the negatives and positives to achieve mass conservation. So cutting the negative values off is not a good idea.

You want to limit sharp horizontal gradients. What is happening at the first rho point next to teh river source? is all that sed being resuspended?

-john

[c.drinkorn]

Thanks for your promt answer!

Speaking of sharp horizontal gradients: I do indeed discover erroneous negative rho values (in fact of the same order as the susp. load concentration) next to the river discharge, so maybe this isn't exactly a sediment issue. T and S look fine, though. The errors occur after several months of simulation and happen to evolve from just one spot. Also, they seem to start rather from the bottom than the surface. Here is my include so you can see which flags I chose. Maybe some of them are prone to instabilities in combination... Unfortunately, I don't have enough experience to spot issues right away but I try hard to understand a little more each day.

Code: Select all

/*
** svn $Id: basin.h 8 2007-02-06 19:00:29Z arango $
*******************************************************************************
** Copyright (c) 2002-2007 The ROMS/TOMS Group
**
**   Licensed under a MIT/X style license
**
**   See License_ROMS.txt
**
*******************************************************************************
**
**  Options for ROMS
*/

#define PROFILE

#define RADIATION_2D 

#undef  DEFLATE 
#define HDF5
#undef PARALLEL_IO
#undef  RST_SINGLE         /* define if single precision restart fields */
#undef  PERFECT_RESTART    /* use to include perfect restart variables */
#define CURVGRID           /* define if using  curvilinear coordinate grid*/
#define GLOBAL_PERIODIC    /* CD: Kate Hedström's patch for global/polar application*/
#undef  DIAGNOSTICS_UV    /* define if writing out dynamics to diagnostics file*/
#undef  DIAGNOSTICS_TS    /* define if writing out tracers  to diagnostics file*/

#define UV_ADV             /* turn ON or OFF advection terms */
#define UV_COR             /* turn ON or OFF Coriolis term */
#define UV_VIS2            /* turn ON or OFF Laplacian horizontal mixing */
#undef  UV_VIS4            /* turn ON or OFF biharmonic horizontal mixing */
#undef  UV_U3ADV_SPLIT     /* use 3rd-order upstream split momentum advection */
#define UV_U3HADVECTION    /* define if 3rd-order upstream horiz. advection */
#undef  UV_SADVECTION      /* turn ON or OFF splines vertical advection */
#undef  UV_C4HADVECTION    /* define if 4th-order centered horizontal advection */
#define UV_QDRAG            /* turn ON or OFF quadratic bottom friction */
#undef UV_LDRAG
#define UV_SMAGORINSKY

#define  VISC_GRID          /* viscosity coefficient scaled by grid size */
#define NONLIN_EOS         /* define if using nonlinear equation of state */
#undef  WJ_GRADP           /* Weighted density Jacobian (Song, 1998) */
#define DJ_GRADPS          /* Splines density  Jacobian (Shchepetkin, 2000) */
#undef  PJ_GRADPQ4         /* quartic 4 Pressure Jacobian (Shchepetkin,2000) */
#define  DIFF_GRID          /* diffusion coefficient scaled by grid size */
#undef  ANA_SPONGE         /* visc_factor / diff_factor from analytical*/

#define TS_DIF2            /* turn ON or OFF Laplacian horizontal mixing */
#undef  TS_DIF4            /* turn ON or OFF biharmonic horizontal mixing */
#undef  TS_U3ADV_SPLIT     /* use 3rd-order upstream split tracer advection */
#define TS_U3HADVECTION    /* define if 3rd-order upstream horiz. advection */
#undef  TS_A4HADVECTION    /* define if 4th-order Akima horiz. advection */
#undef  TS_C4HADVECTION    /* define if 4th-order centered horizontal advection */

#undef  TS_MPDATA          /* define if recursive MPDATA 3D advection */

#define TS_A4VADVECTION    /* define if 4th-order Akima vertical advection */
#undef  TS_C4VADVECTION    /* define if 4th-order centered vertical advection */
#undef  TS_SVADVECTION     /* define if splines vertical advection */
#undef  TS_SMAGORINSKY     /* define if Smagorinsky-like diffusion */

#undef  MIX_S_TS           /* mixing on constant S-surfaces */
#define MIX_GEO_TS         /* mixing on geopotential (constant Z) surfaces */
#define MIX_S_UV           /* mixing along constant S-surfaces */
#undef  MIX_GEO_UV         /* mixing on geopotential (constant Z) surfaces */
#define SALINITY           /* define if using salinity */
#define SOLVE3D            /* define if solving 3D primitive equations */
#undef  BODYFORCE          /* define if applying stresses as bodyforces */
#undef  SPLINES            /* turn ON or OFF parabolic splines reconstruction */
#define MASKING            /* define if there is land in the domain */
#define AVERAGES           /* define if writing out time-averaged data */
#undef  AVERAGES_DETIDE    /*use if writing out NLM time-averaged detided fields*/
#undef  AVERAGES_AKT       /*undef Frode, define if writing out time-averaged AKt*/
#undef  AVERAGES_AKS       /*undef Frode, define if writing out time-averaged AKs*/
#undef  AVERAGES_AKV       /* define if writing out time-averaged AKv */

#undef  STATIONS           /* define if writing out station data */
#undef  STATIONS_CGRID     /* define if extracting data at native C-grid */

#undef  BVF_MIXING         /* define if Brunt_Vaisala frequency mixing */
#undef  LMD_MIXING         /* define if Large et al. (1994) interior closure */
#undef  MY25_MIXING        /* define if Mellor/Yamada level-2.5 mixing */
#define  GLS_MIXING         /* Activate Generic Length-Scale mixing */

#ifdef  GLS_MIXING
# define N2S2_HORAVG       /* Activate horizontal smoothing of buoyancy/shear */
# undef  KANTHA_CLAYSON    /* Value for CLS_CMU0 and CLS_C3M vary with choise of stability function */
# define CANUTO_A           
# undef  CANUTO_B
#endif
#ifdef LMD_MIXING
# define  LMD_BKPP          /* use if bottom boundary layer KPP mixing */
# define  LMD_CONVEC        /* use to add convective mixing due to shear instability */
# define  LMD_DDMIX         /* use to add double-diffusive mixing */
# define  LMD_NONLOCAL      /* use if nonlocal transport */
# define  LMD_RIMIX         /* use to add diffusivity due to shear instability */
# define  LMD_SHAPIRO       /* use if Shapiro filtering boundary layer depth */
# define  LMD_SKPP          /* use if surface boundary layer KPP mixing */
#endif

#define ANA_BSFLUX         /* analytical bottom salinity flux */
#define ANA_BTFLUX         /* analytical bottom temperature flux */
#undef  ANA_GRID           /* analytical grid set-up */
#undef  ANA_INITIAL        /* analytical initial conditions */
#undef  ANA_MASK           /* analytical land/sea masking */
#undef  ANA_MEANRHO         
#undef  ANA_SSFLUX         /* analytical surface salinity flux */
#undef  ANA_STFLUX         /* analytical surface temperature flux */
#undef  ANA_SMFLUX         /* analytical surface momentum stress */
#undef  ANA_NUDGCOEF

/* CLIMATOLOGY */
#undef  M2CLIMATOLOGY      /* define 2D momentum climatology */
#undef  M3CLIMATOLOGY      /* define 3D momentum climatology */
#undef  TCLIMATOLOGY       /* define tracers climatology */
#undef  M2CLM_NUDGING
#undef  M3CLM_NUDGING      /* nudging 3D momentum climatology */
#undef  TCLM_NUDGING       /* nudging tracers climatology */
#undef  OBC_NUDGING

#undef  WRF_COUPLING       /* coupling to WRF atmospheric model */

/* ATMOSPHERIC FORCING */
#define BULK_FLUXES        /* turn ON or OFF bulk fluxes computation */
#ifdef  BULK_FLUXES
# undef  ANA_RAIN          /* analytical rain fall rate */
# undef  ANA_PAIR          /* analytical surface air pressure */
# undef  ANA_HUMIDITY      /* analytical surface air humidity */
# undef  ANA_CLOUD         /* analytical cloud fraction */
# undef  ANA_TAIR          /* analytical surface air temperature */
# undef  ANA_WINDS         /* analytical surface winds */
# define EMINUSP           /* turn ON internal calculation of E-P */
# define ANA_SRFLUX        /* analytical surface shortwave radiation flux */
# define ALBEDO            /* use albedo equation for shortwave radiation */
# undef  CLOUDS            /* This is obsolete since globaldefs switches it on again*/
# undef  LONGWAVE_OUT      /* compute outgoing longwave radiation */
# define LONGWAVE          /* Compute net longwave radiation internally */
# define COOL_SKIN         /* turn ON or OFF cool skin correction */
# undef  SHORTWAVE
#endif
#define ATM_PRESS          /* use to impose atmospheric pressure onto sea surface */
#define SOLAR_SOURCE       /* define solar radiation source term */
#undef SPECIFIC_HUMIDITY  /* if input is specific humidity in kg/kg */

/* SEDIMENTATION */
#define SEDIMENT           /* switch on sedimentation module */
#define ANA_SEDIMENT       /* analytical initial properties of sediments */
#define SED_MORPH          /* evolving model elevation */
#define SUSPLOAD           /* Suspended load transport */
#define SED_DENS           /* Let sediment affect equation of state */
#define ANA_SPFLUX         /* Surface passive tracer flux from analytical */
#define ANA_BPFLUX         /* Bottom passive tracer flux from analytical */

/* TIDES */
#undef SSH_TIDES          /* turn on computation of tidal elevation, default define */
#undef UV_TIDES           /* turn on computation of tidal currents */
#undef ADD_FSOBC          /* Add tidal elevation to processed OBC data, default define */
#undef ADD_M2OBC          /* Add tidal currents  to processed OBC data */
#undef ADD_FS_INV_BARO
#undef RAMP_TIDES         /* Spin up tidal forcing */


/* ELVER */
#undef UV_PSOURCE         /* turn ON or OFF point Sources/Sinks */
#undef TS_PSOURCE         /* turn ON or OFF point Sources/Sinks */
/* --------------------------- */

[jcwarner]

so a suggestion would be to change
define TS_U3HADVECTION
to
define TS_MPDATA
and see how that looks.

Also, just a question: are the river sources all putting flow into the domain? some people use tidal river forcing that oscillates +/- to simulate a tidal flow.
when you say negative sediment next to the river discharge. is that the first rho point adjacent to the source? does bed_mass go to 0 before the negative conc? trying to think of some ways to help.

[c.drinkorn]

MPDATA did really solve the problem which does of course make sense because it just restricts the local simulation results. Even though I am always a bit concerned with too much correction terms in the calculation, I think switching to MPDATA was a really good idea. Thank you for this great advice!

Regarding your other questions: I have only momentum and tracer (T and S, no discharge - I am still considering this, though) transport from the rivers. No volume (ergo water mass) is being introduced.
The negative concentrations seemingly didn't appear right at the first rho point from the river point source but rather in their vicinity. And I don't extract bed mass so far (reducing simulation time) but bed thickness was not much affected at the areas of question...
So I guess it was more some kind of instability that the MPDATA scheme is capable of constraining. It would have taken me forever to find this myself.

Thank you again!

[c.drinkorn]

Oh and by the way, I also needed to undefine the vertical tracer advection scheme, of course. Just in case someone is reading this later...

[c.drinkorn]

Something else: Is there anything I need to switch on in order to write out bed mass? Setting the true flag in sediment.in didn't do the trick.

[jcwarner]

i get bed mass written out.
make sure you set it for sand or mud, whichever you have activated.

[c.drinkorn]

Mysterious. I think, I have everything set as it should:

Code: Select all

...
Hout(idmud)  == T       ! mud_01, ...             suspended concentration
Hout(iMfrac) == T       ! mudfrac_01, ...         bed layer fraction
Hout(iMmass) == T       ! mudmass_01, ...         bed layer mass
Hout(iMUbld) == F       ! bedload_Umud_01, ...    bed load at U-points
Hout(iMVbld) == F       ! bedload_Vmud_01, ...    bed load at V-points
...

I also tried to switch on bed layer fraction in order to test this one and it gets written out.

[jcwarner]

are you doing mud or sand?
NNS =
NCS =

[c.drinkorn]

MUd, so NCS equals 1. I thought I found the reason but still no writing of bed mass:
In the Write-out section of sediment_inp.h it says

Code: Select all

...
              DO itrc=1,NST
                i=idfrac(itrc)
                IF (Hout(i,ng)) WRITE (out,160) Hout(i,ng),             &
     &              'Hout(idfrac)',                                     &
     &              'Write out bed fraction, sediment ', itrc,          &
     &              TRIM(Vname(1,i))
              END DO
              DO itrc=1,NST
                i=idBmas(itrc)
                IF (Hout(i,ng)) WRITE (out,160) Hout(i,ng),             &
     &              'Hout(idfrac)',                                     &
     &              'Write out mass, sediment ', itrc,                  &
     &              TRIM(Vname(1,i))
              END DO
...

I thought it must be "idBmas" instead of "idfrac" for the MASS instead. But this didn't solve the problem. In the output there is no mention of writing out bed mass. So somehow it is simply not processed and therefore doesn't matter what it says in the bed mass case.

[jcwarner]

does it report in the std out at the top with all the other zob, tcline, etc etc stuff, then it gets to sediment, it should report the sed classes, then does it say it will Hout().... and write out the sediment mass??

[c.drinkorn]

It does define mud as sediment class but not mention bed mass as Hout output:

Code: Select all

 Sediment Parameters, Grid: 01
 =============================


 Size     Sd50        Csed        Srho        Wsed        Erate       poros
 Class    (mm)       (kg/m3)     (kg/m3)     (mm/s)     (kg/m2/s)    (nondim)

   1    0.0000E+00  0.0000E+00  1.8000E+03  1.0000E-01  5.0000E-04  9.0000E-01

         tau_ce      tau_cd      nl_tnu2     nl_tnu4     Akt_bak      Tnudg
         (N/m2)      (N/m2)      (m2/s)      (m4/s)       (m2/s)      (day)

   1    1.0000E-01  1.0000E-01  0.0000E+00  0.0000E+00  5.0000E-06  0.0000E+00

         morph_fac
         (nondim)

   1    1.0000E+00

 New bed layer formed when deposition exceeds  0.10000E-01 (m).
 Two first layers are combined when 2nd layer smaller than  0.00000E+00 (m).
 Rate coefficient for bed load transport =  0.50000E-01

          F  LtracerSponge(03Turning OFF sponge on tracer 03: mud_01
          F  LtracerSrc(03)  Turning OFF point sources/Sink on tracer 03: mud_01
          F  LtracerCLM(03)  Turning OFF processing of climatology tracer 03: mud_01
          F  LnudgeTCLM(03)  Turning OFF nudging of climatology tracer 03: mud_01

          T  Hout(idTvar)   Write out sediment01: mud_01
          T  Hout(idfrac)   Write out bed fraction, sediment 01: mudfrac_01
          T  Hout(idSbed)   Write out BED property 01: bed_thickness
          T  Hout(idSbed)   Write out BED property 02: bed_age

          T  Aout(idTvar)   Write out averaged sediment01: mud_01

I looked into all routines that could possibly be processing bed mass but I found nothing (except for the one in sediment_inp.h). I am lost here.

[jcwarner]

can you send your sediment.in

[c.drinkorn]

Ok, here is my sediment.in

[c.drinkorn]

In my restart files I get mudmass written out so in principle this works. Just not in the history file...

[jcwarner]

think i found it!

in sediment_inp.h there is a section of code
CASE ('Qout(iMmass)')
Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
DO ng=1,Ngrids
DO itrc=1,NCS
i=idBmas(itrc)
Hout(i,ng)=Lmud(itrc,ng)
END DO
END DO

but that 6th line should be
Qout(i,ng)=Lmud(itrc,ng)

So the Hout was being overwritten. Can you find that section of code and modify
Hout(i,ng)=Lmud(itrc,ng)
to
Qout(i,ng)=Lmud(itrc,ng)

-j

[c.drinkorn]

Yessss, that was it! I didn't look in the Qout section. Thank you!