12.2.1. USGS Sediment Model#

This USGS sediment model is derived from the UCLA/USGS ROMS community. See Blaas et al. [2007], Warner et al. [2008] and Shafiei [2021] for details.

Regarding the time and space resolution considered, the explicit solution generally refers to quantities averaged over wave periods, although the implementation of a nonhydrostatic solver in CROCO opens the way to a wave-resolved approach. One of the crucial ingredients in the sediment transport model is a reliable representation of wave-averaged (or wave-resolved) hydrodynamics and turbulence.

In the wave-averaged approach, the wave boundary layer is not resolved explicitly, but the lower part of the velocity and sediment concentration profile in the current boundary layer is important for the calculation of the sediment transport rates. Similarly, an accurate assessment of the bottom boundary shear stress (including wave effects) is required since it determines the initiation of grain motion and settling and resuspension of suspended load (see BBL). Thus, the sediment concentration and current velocity profiles in the unresolved part of the near-bottom layer have to be parameterized. Characterization of the sediments (mainly density and grain size, making general assumptions about shape and cohesiveness) is done either as a time-dependent prescribed function at the point sources or at the sea bed as an initial (soon space-dependent) condition. Sediment concentration may be considered as passive with respect to the flow density or as active if concentration values require such (the latter is not implemented yet).

12.2.1.1. Sediment bed#

The sediment bed is represented by three-dimensional arrays with a fixed number of layers beneath each horizontal model cell. Each cell of each layer in the bed is initialized with a thickness, sediment-class distribution, porosity, and age. The mass of each sediment class in each cell can be determined from these values and the grain density. The bed framework also includes two-dimensional arrays that describe the evolving properties of the seabed, including bulk properties of the surface layer (active layer thickness, mean grain diameter, mean density, mean settling velocity, mean critical stress for erosion) and descriptions of the subgrid scale morphology (ripple height and wavelength). These properties are used to estimate bed roughness in the BBL formulations and feed into the bottom stress calculations. The bottom stresses are then used by the sediment routines to determine resuspension and transport, providing a feedback from the sediment dynamics to the hydrodynamics.

The bed layers are modified at each time step to account for erosion and deposition and track stratigraphy. At the beginning of each time step, an active layer thickness \(z_a\) is calculated [Harris and Wiberg, 1997]. \(z_a\) is the minimum thickness of the top bed layer. If the top layer is thicker than \(z_a\), no action is required. If the top layer is less than \(z_a\), then the top layer thickness is increased by entraining sediment mass from deeper layers until the top layer thickness equals \(z_a\). If sediment from deeper than the second layer is mixed into the top layer, the bottom layer is split to enforce a constant number of layers and conservation of sediment mass. Each sediment class can be transported by suspended-load and/or bedload (below). Suspended-load mass is exchanged vertically between the water column and the top bed layer. Mass of each sediment class available for transport is limited to the mass available in the active layer. Bedload mass is exchanged horizontally within the top layer of the bed. Mass of each sediment class available for transport is limited to the mass available in the top layer. Suspended-sediment that is deposited, or bedload that is transported into a computational cell, is added to the top bed layer. If continuous deposition results in a top layer thicker than a user-defined threshold, a new layer is provided to begin accumulation of depositing mass. The bottom two layers are then combined to conserve the number of layers. After erosion and deposition have been calculated, the active-layer thickness is recalculated and bed layers readjusted to accommodate it. This step mixes away any very thin layer (less than the active layer thickness) of newly deposited material. Finally the surficial sediment characteristics, such as D50, ripple geometry, etc., are updated and made available to the bottom stress calculations.

12.2.1.2. Suspended-sediment transport#

The concentration of sediment suspended in the water column is transported, like other conservative tracers (e.g., temperature and salinity) by solving the advection–diffusion equation with a source/sink term for vertical settling and erosion:

\[ \underbrace{\frac{\partial C }{\partial t} }_{RATE} = - \underbrace{\vec{\nabla} . \vec{\textbf v} C }_{ADVECTION} + \underbrace{\mathcal{D}_C }_{MIXING} - \underbrace{\frac{\partial w_s C }{\partial z} }_{SETTLING} + \underbrace{ \frac{ E }{\delta z_b}\Big|_{z=z_b} }_{EROSION}\]

\(C\) is the Reynolds-averaged, wave-averaged (unless used in wave-resolving mode) sediment concentration of a particular size class; \(\vec{\textbf v}\) is the flow velocity (it is the Lagrangian velocity \(\vec{\textbf v_L}\) in wave-averaged equations, comprising the Stokes drift \(\vec{\textbf v_S}\)).

For each size class, the source or sink term represents the net of upward flux of eroded material E and downward settling, i.e., the deposition flux. \(w_s\) is the settling velocity, dependent on sediment grain size, but independent of flow conditions and concentrations. It is an input parameter of the model (WSED in sediment.in; see below). Settling is computed via a semi-Lagrangian advective flux algorithm, which is unconditionally stable [Durran, 2010]. It uses a piece-wise parabolic vertical reconstruction of the suspended sediment for high-order interpolation, with WENO constraints to avoid oscillations. \(E\) is the erosion flux at the sea floor and is only applied to the first grid level of height \(z_b\) and cell size \(\delta z_b\). The erosion flux for each class is given by:

\[E = E_0 (1-p) \, \phi \left( \frac{\tau_s}{\tau_{c}} -1 \right) ~\text{for} ~\tau_s > \tau_{c}\]

\(E_0\) is an empirical erosion rate (ERATE parameter in sediment.in; see below); p is the sediment porosity; \(\phi\) is the volumetric fraction of sediment of the class considered; \(\tau_{c}\) is the critical shear stress; and \(\tau_s\) is the shear stress magnitude on the grains (skin stress due to wave-induced bed orbital velocities and mean bottom currents; see BBL). The critical shear stress is the threshold for the initiation of sediment motion.

Zero-flux boundary conditions are imposed at the surface and bottom in the vertical diffusion equation. Lateral open boundaries are treated as other tracers according to Marchesiello et al. [2001]. A quasi-monotonic 5th-order advection scheme (WENO5-Z, Borges et al. [2008]) can be used for horizontal and vertical advection of all tracers, including sediments.

12.2.1.3. Bedload transport#

The bedload flux \(q_b\), which is considered unresolved by the model can be calculated using different bedload models implemented in CROCO. The formulation by Meyer-Peter Muller [Meyer-Peter and Müller, 1948] is suited to rivers or continental shelf problems, where nonlinear wave effects are small. For nearshore applications, where wave nonlinearity is important, the bedload transport formulation proposed by van der A et al. (2013) is implemented as in Shafiei [2021] following Kalra et al. [2019] with some modifications.

Each formulation depends on the characteristics of individual sediment classes, including median size \(d_{50}\), grain density \(\rho_s\), specific density in water \(s=\rho/\rho_s\), and critical shear stress \(\tau_{c}\). Non-dimensional transport rates \(\Phi\) are calculated for each sediment class and converted to dimensional bedload transport rates \(q_b\) using:

\[q_b = \Phi \sqrt{ (s-1) g d_{50}^3 \rho_s }\]

These are horizontal vector quantities with directions that correspond to the combined bed-stress vectors. Details on the computation of \(\Phi\) differs in the Meyer-Peter Müeller or van der A formulations.

Slope effect: bedload fluxes are corrected to account for the avalanche process, i.e., the gravitational flow of sand occuring when the bottom slope exceeds the critical slope angle:

\[q_{b,slope} = q_{b} \left( \dfrac{0.65} {(0.65-\tan\beta)\cos\beta}\right),\]

This correction considers the effect of the bed slope \(\beta=\tan^{-1}(d z_b / dx)\). The value 0.65 is derived from the consideration of an angle of repose of \(33^\circ\).

Bedload numerics: bedload fluxes are computed at grid-cell centers and are limited by the availability of each sediment class in the top layer. Fluxes are then interpolated on cell faces using an upwind approach, either 1rst-order (e.g., Lesser et al. [2004]) or 5th order, or even a WENO5 interpolation to avoid oscillations. Flux differences are then used to determine changes of sediment mass in the bed at each grid cell.

12.2.1.3.1. Meyer-Peter Müller : Transport by currents#

Meyer-Peter Müller [Meyer-Peter and Müller, 1948] formulation :

\[\Phi= max \left [ 8 (\theta_s-\theta_c)^{1.5},0 \right ]\]

where \(\Phi\) is the magnitude of the non-dimensional transport rate for each sediment class, \(\theta_s\) is the non-dimensional Shields parameter for skin stress:

\[\theta_s=\frac{\tau_s}{(s-1)g d_{50}}\]

\(\theta_c\) is the critical Shields parameter, and \(\tau_s\) the magnitude of total skin-friction component of bottom stress computed from:

\[\tau_s=\sqrt{\tau_{sx}^2 + \tau_{sy}^2}\]

where \(\tau_{sx}\) and \(\tau_{sy}\) are the skin-friction components of bed stress, from currents alone or the maximum wave-current combined stress, in the x and y directions. These are computed at cell faces (u and v locations) and then interpolated to cell centers ( \(\rho\) points). The bedload transport vectors are partitioned into x and y components based on the magnitude of the bed shear stress as:

\[\begin{split}q_{bx}=q_{b}\frac{\tau_{sx}}{\tau_{s}} \\ q_{by}=q_{b}\frac{\tau_{sy}}{\tau_{s}}\end{split}\]

12.2.1.3.2. van der A (2013): Transport by nonlinear waves#

The SANTOSS bedload model is based on the half-wave cycle concept proposed by Dibajnia and Watanabe [1993] that captures asymmetric transport by non-linear waves and the effect of phase lag between mobilization and transport. CROCO contains an adapted version by Shafiei [2021], based on the implementation of Kalra et al. [2019]. In our formulation, the effect of wave-averaged currents is removed by default (assuming that the transport by currents is effectively performed by the suspended load model) and it thus only retains the nonlinear effects of waves. In the brief presentation below, we retain the current terms for completeness, but focus on wave effects.

The method to obtain bedload transport under asymmetric waves can be divided into three major steps (van der A, 2013) [Kalra et al., 2019]. In the first step, the asymmetric waveform based on the Ursell number is evaluated using wave statistics. The Shields parameter for each half cycle of the wave form is computed in the second step. Finally, a phase lag is estimated from the velocity and sediment concentrations that determine the amount of bedload transported in the half cycle following mobilization. The non-dimensional bedload transport rate \(\Phi\) is thus given by:

\[\Phi = \dfrac{1}{T}\left[ \dfrac{\theta_{c}}{\big\vert\theta_{c}\big\vert^{1/2}} T_c \left( \Omega_{cc}+\dfrac{T_c}{2T_{cu}} \Omega_{tc}\right) + \dfrac{\theta_{t}}{\big\vert\theta_{t}\big\vert^{1/2}} T_t \left( \Omega_{tt}+\dfrac{T_t}{2T_{tu}} \Omega_{ct}\right) \right],\]

where \(T, T_c, T_t, T_{cu}\) and \(T_{tu}\) are the wave period, duration of wave crest half cycle, duration of wave trough half cycle, duration of accelerating flow within the crest half cycle and duration of accelerating flow within the trough half cycle respectively, \(\theta_{c}\) and \(\theta_{t}\) represent the Shields numbers associated with the wave crest and trough half cycles. The sand load transported during the crest period is the combination of \(\Omega_{cc}\) (mobilized during the crest period) and \(\Omega_{tc}\) (mobilized during the trough period). Similarly, \(\Omega_{tt}\) and \(\Omega_{ct}\) are the sand load transported during the trough period (mobilized during the trough and crest periods respectively).

The sand load transported during each half-cycle is conventionally modeled according to a power law of Shields number:

\[\Omega_{i} = \max\left( 11 \left(\big\vert\theta_{i}\big\vert - \theta_{cr}\right)^{1.2},0\right),\]

where \(\theta_{cr}\) is the critical Shields number and, hereafter, the subscript “i” is either “c” for crest or “t” for trough half cycles. To determine \(\Omega_{ct}\) and \(\Omega_{tc}\), i.e., the portion of the bedload remaining in suspension to be transported in the next half cycle, a phase lag parameter is evaluated.

Let’s assume a two-dimensional (x,z) cross-shore problem for simplicity. The Shields number for the peak or trough (\(\theta_{i}=\theta_{t}\) or \(\theta_{c}\)) is calculated according to:

\[\theta_{i} = \dfrac{ \tfrac{1}{2} f_{w\delta i} \vert u_{i,r}\vert u_{i,r} }{(s-1) g d_{50}}.\]

\(u_{i,r}\) is the representative cross-shore combined wave-current velocity at trough or crest half cycles calculated as:

\[u_{i,r}=\dfrac{\hat u_{i}}{\sqrt 2} + \vert u_{\delta} \vert,\]

where \(\hat u_{i}\) is the peak crest or trough orbital velocities, \(u_{\delta}\) is the steady current velocity at the top of the wave boundary layer. \(f_{w\delta i}\) is the linear wave-current friction factor at crest or trough calculated by Ribberink [1998]:

\[f_{w\delta i}= \dfrac{\hat{u}}{u_\delta + \hat{u}} f_{w i} + \dfrac{u_\delta}{u_\delta + \hat{u}} f_\delta\]

where \(\hat{u}\) is the representative orbital velocity amplitude for the whole flow cycle (given by \(\hat{u}=\sqrt{2} u_{orb}\)). \(f_\delta\) is the current-related friction factor dependent on a current-related roughness \(k_{s\delta}\) and \(f_{wi}\) is the wave friction factor, calculated separately for the crest and trough half-cycles and depends on a wave-related roughness \(k_{sw}\). If the representative orbital excursion amplitude \(\hat{a}=\hat{u}T / 2\pi\) is large enough (i.e., greater than 1.587 \(k_{sw}\)):

\[f_{wi} = 0.00251 \, e^{5.21 \left[ \left(\tfrac{2 T_{iu}}{T_i} \right)^{2.6} \, \tfrac{\hat{a}}{k_{sw}} \right]^{-0.19}},\]

otherwise, \(f_{wi}=0.3\).

If \(u_\delta=0\) (\(u_{i,r}=\hat u_{i}/\sqrt 2\) and \(f_{w\delta i}=f_{wi}\)), the effect of currents is completely removed from the bedload transport calculation. This choice represents our default to avoid double counting the transport by wave-averaged currents.

Finally, following Kalra et al. [2019], after calculating \(\Phi\), we apply to \(q_{b}\) a bedload factor \(f_{bld}\) (bedload_coeff in CROCO). This factor allows us to adjust the relative contribution to sediment transport of wave-induced bedload compared to the suspended load transported by mean currents. In this way, we have a better control of the antagonistic mechanisms that govern onshore and offshore transports respectively.

12.2.1.4. Morphology#

The bed evolution (variation in time of \(z_b\), the height of the bed), is calculated from the divergence of sediment fluxes (Exner equation), which results from the difference between erosion and sedimentation of suspended sediments. In wave-averaged equations, where residual wave effects need to be parametrized as bedload fluxes \(q_b\), the bed evolution also arises from the divergence of these fluxes.

\[\dfrac{\partial z_b}{\partial t} = - \dfrac{f_{mor}}{1-p} \left( \dfrac{\partial q_b} {\partial x} - w_s \dfrac{\partial C}{\partial z} + E \right).\]

This equation accounts for a morphological acceleration factor \(f_{mor}\) (morph_fac in CROCO). A value of 1 has no effect, and values greater than 1 accelerate the bed response. The concept of morphological acceleration is based on the fact that morphodynamic changes are slower than hydrodynamic ones [van Rijn, 1993]. In this case, the bed evolution can be accelerated without affecting the hydro-morphological solution. The increased rate of morphological change can be useful for simulating evolution over long time periods. Strategies for morphological updating are described by Roelvink [2006] and implemented in CROCO following Warner et al. [2008]. In our implementation, bedload fluxes, erosion, and deposition rates are multiplied by \(f_{mor}\), while the magnitude of sediment concentrations in the water column is not modified – just the exchange rate to and from the bed. For both bedload and suspended load, sediment is limited in availability, based on the true amount of sediment mass (not multiplied by the scale factor).

For dynamical consistency, the vertical velocity is modified (in omega.F) by the rate of change of vertical grid levels \(dz/dt\), adjusting to the moving sea floor and free surface (grid “breathing” component; Shchepetkin and McWilliams [2005]). This method is mass conserving and retains tracer constancy preservation.

12.2.1.5. Sediment Density#

Witt CPP SED_DENS, effects of suspended sediment on the density field are included with terms for the weight of each sediment class in the equation of state for seawater density as:

\[\rho = \rho_{water}+\sum_{m=1}^{Nsed}\frac{C_m}{\rho_{s,m}}\left(\rho_{s,m}-\rho_{water}\right)\]

This enables the model to simulate processes where sediment density influences hydrodynamics, such as density stratification and gravitationally driven flows.

Related CPP options:

SUSPLOAD	Activate suspended load transport
BEDLOAD	Activate bedload transport
MORPHODYN	Activate morphodynamics
BEDLOAD_VANDERA	van der A formulation for bedload (van der A et al., 2013)
BEDLOAD_MPM	Meyer-Peter-Muller formulation for bedload [Meyer-Peter and Müller, 1948]
SLOPE_LESSER	Lesser formulation for avalanching [Lesser et al., 2004]
SLOPE_NEMETH	Nemeth formulation for avalanching (Nemeth et al, 2006)
BEDLOAD_UP1	Bedload flux interpolation: upwind 1rst order
BEDLOAD_UP5	Bedload flux interpolation: upwind 5th order
BEDLOAD_WENO5	Bedload flux interpolation: WENO 5th order
ANA_SEDIMENT	Set analytical sediment size, initial ripple and bed parameters
ANA_BPFLUX	Set kinematic bottom flux of sediment tracer (if different from 0)
SPONGE_SED	Gradually reduce erosion/deposition near open boundaries
SED_DENS	Activate the effect of suspended sediment on the density field

Preselected options:

#ifdef SEDIMENT
# undef  MUSTANG
# define ANA_SEDIMENT
# define SPONGE_SED
# define Z0_BL
# define Z0_RIP
# ifdef BEDLOAD
#  ifdef BEDLOAD_VANDERA       /* default BEDLOAD scheme  */
#  elif defined BEDLOAD_MPM
#  elif defined BEDLOAD_WULIN
#  elif defined BEDLOAD_MARIEU
#  else
#   if (defined WAVE_OFFLINE || defined WKB_WWAVE ||\
        defined ANA_WWAVE    || defined OW_COUPLING)
#    define BEDLOAD_VANDERA
#   else
#    define BEDLOAD_MPM
#   endif
#  endif
#  ifdef BEDLOAD_UP1           /* default INTERPOLATION  */
#  elif defined BEDLOAD_UP5
#  elif defined BEDLOAD_WENO5
#  else
#   define BEDLOAD_UP1
#  endif
#  ifdef SLOPE_LESSER          /* default SLOPE scheme */
#  elif defined SLOPE_NEMETH
#  elif defined SLOPE_KIRWAN
#  else
#   define SLOPE_LESSER
#  endif
# endif /* BEDLOAD */
 #endif /* SEDIMENT */

Parameters in sediment.in

Stitle (a80)
CROCO - Sediment - Test

Sd(1-NST), CSED, SRHO,  WSED,  ERATE,  TAU_CE, TAU_CD, BED_FRAC(1:NLAY)
125    9.9   2650.  9.4    25.0e-5  0.05    0.14    0.4  0.4
050    0.0   2650.  1.6     4.0e-5  0.01    0.14    0.6  0.6

BTHK(1:NLAY)
   10.

BPOR(1:NLAY)
41  0.42

Hrip
03

Lrip
14

bedload_coeff


morph_fac


transC
03

transN
2

tcr_min
03

tcr_max
5

tcr_slp
3

tcr_off


tcr_tim


L_ADS  L_ASH  L_COLLFRAG  L_TESTCASE
    F       T         F           F

F_DP0      F_NF    F_DMAX     F_NB_FRAG    F_ALPHA   F_BETA  F_ATER   F_ERO_FAC F_ERO_NBFRAG  F_COLLFRAGPARAM  F_CLIM    F_ERO_IV
000004    2.     0.0015        2.          0.35     0.15    0.         0.        2.             0.01          0.001        1

MUD_FRAC_EQ [1:NMUD]
10  0.20  0.40  0.20  0.10  0.0  0.0  0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

MUD_T_DFLOC


END of sediment input data

GLOSSARY

CARD 1: String with a maximum of eighty characters.

Stitle : Sediment case title.

CARD 2: Sediment grain parameters & initial values (NST lines)

Sd : Diameter of grain size class [mm].
CSED : Initial concentration (spatially uniform) [kg/m3].
SRHO : Density of sediment material of size class [kg/m3]. Quartz: SRHO=2650 \(~kg/m^3\)
WSED : Settling velocity of size class [mm/s].

Typically [Soulsby, 1997]:

WSED = \(10^3 ~(visc ~(\sqrt{10.36^2 + 1.049 D^3}~-10.36)~/~D_{50}\) [mm/s]

with
- \(D=D_{50} ~(g ~(\mathrm{SRHO}/\rho_0-1)~/~(visc^2)~)^{0.33333}\)
- \(D_{50}=10^{-3} Sd\) [m]
- \(visc=1.3 ~10^{-3}/\rho_0\) [m2/s]
ERATE : Erosion rate of size class [kg/m2/s].

Typically:

\(\text{ERATE}=10^{-3} ~\gamma_0 ~\mathrm{WSED} ~\mathrm{SRHO} ~~[kg/m^2/s]\)

with \(\gamma_0= 10^{-3}-10^{-5}\) [Smith and McLean, 1977]
TAU_CE : Critical shear stress for sediment motion [N/m2] (initiation of bedload for coarses, suspension for fines).

Typically :

\(\mathrm{TAU_CE}= 6.4~10^{-7} ~\rho_0 ~\mathrm{WSED}^2\) [N/m2]
TAU_CD : Critical shear stress for deposition of cohesive sediments [N/m2]
BED_FRAC : Volume fraction of each size class in each bed layer (NLAY columns)

[0<BED_FRAC<1]

CARD 3: Sediment bed thickness, 1st field is top layer (‘delt_a’)

BTHK : Initial thicknesses of bed layers [m] Bthk(1) active layer thickness, fixed in simulation unless SUM(Bthk(:))<Bthk(1)

CARD 4: Sediment bed porosity

BPOR : Initial porosity of bed layers [m] used in ana_sediment ifdef ANA_SEDIMENT (not in init.nc)

CARD 5: Bottom ripple height

Hrip : Initial ripple height [m] used in ana_sediment ifdef ANA_SEDIMENT (not in init.nc)

CARD 6: Bottom ripple length

Lrip : Initial ripple length [m] used in ana_sediment ifdef ANA_SEDIMENT (not in init.nc)

CARD 7: Bedload coefficient

bedload_coeff : factor limiting the magnitude of bedload flux 0<bedload_coef<1

CARD 8: Morphological acceleration factor

morph_fac : factor accelerating bed evolution morph_fac>=1

CARD 9 :

transC : Cohesive transition- Under that value of total mud fraction entire bed behaves as a non-cohesive bed

CARD 10 :

transN : Noncohesive transition- Over that value of total mud fraction entire bed behaves as a cohesive bed

CARD 11 :

tcr_min : Minimum shear for erosion

CARD 12 :

tcr_max : Maximum shear for erosion

CARD 13 :

tcr_slp : Tau_crit profile slope

CARD 14 :

tcr_off : Tau_crit profile offset

CARD 15 :

tcr_tim : Tau_crit consolidation rate

CARD 16 : booleans for flocculation

L_ADS : Boolean set to .true. if differential settling aggregation
L_ASH : Boolean set to .true. if shear aggregation
L_COLLFRAG : Boolean set to .true. if collision-induced fragmentation enable
L_TESTCASE : If .TRUE. sets G(t) to values from Verney et al. [2011] lab experiment

CARD 17 : flocculation Sediment Parameters

F_DP0 : Primary particle size (m), typically 4e-6 m
F_NF : Floc fractal dimension, typically ranging from 1.6 to 2.6
F_DMAX : Maximum diameter (m)

USGS Sediment Model

Contents