Sediment Transport 1

From CIRPwiki
Revision as of 15:32, 6 December 2022 by Rdchllbm (talk | contribs) (Undo revision 13527 by Rdchllbm (talk))
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Sediment Transport Overview

For sand transport, the wash-load (i.e. sediment transport which does not contribute to the bed-material) can be assumed to be zero, and therefore, the total-load transport is equal to the sum of the bed- and the suspended-load transports: qtk*=qsk*+qbk*.

There are currently three sediment transport models available in CMS:

(1) Equilibrium total load

(2) Equilibrium bed load plus non-equilibrium suspended load, and

(3) Non-equilibrium total-load.

The first two models are single-size sediment transport models and are only available with the explicit time-stepping schemes. The third is multiple-sized sediment transport model and is available with both the explicit and implicit time-stepping schemes.

Equilibrium Total-load Transport Model

In this model, both the bed load and suspended load are assumed to be in equilibrium. The bed change is solved using a simple mass balance equation known as the Exner equation.

ρs(1ρm')fmorphzbt=qt*jxj+xj(Dsqt*zbxj) (2-42)

for j=1,2;k=1,2,.....,N , where N is the number of sediment size classes and

t = time [s]


h = total water depth [m]
xj =Cartesian coordinate in the jth direction [m]
qt* = equilibrium total-load transport rate [kg/m/s]
zb = bed elevation with respect to the vertical datum [m]
pm' = bed porosity [-]
fmorph = morphologic acceleration factor [-]
ρs = sediment density [~2650 kg/m3 for quartz sediment]
Ds = empirical bed-slope coefficient (constant) [-]

Because the model assumes that both the sediment transport is equilibri-um, it only recommended for coarse grids with resolutions larger than 50-100 m where the assumption of equilibrium sediment transport is more appropriate. As mentioned above the equilibrium total-load sediment transport model is a single-size sediment transport model and is only available with the explicit time-stepping scheme. For more information on the equilibrium sediment transport model, the reader is referred to Buttolph et al. (2007).

Equilibrium Bed-load plus Nonequilibrium Suspended Load Transport Model

Calculations of suspended load and bed load are conducted separately. The bed load is assumed to be in equilibrium and is included in the bed change equation while the suspended load is solved through the solution of an advection-diffusion equation. Actually the advection diffusion equation is a non-equilibrium formulation, but because the bed load is assumed to be in equilibrium, this model is referred to the "Equilibrium A-D" model.

Suspended-load Transport Equation

The transport equation for the suspended load is given by

(hC)t+(hVjC)xj=xj(vshCxj)+EbDb (2-43)

where

t = time[s]
h = water depth [m]
xj = Cartesian coordinate in the
Eb entrainment or pick-up function [kg/m2/s]
Db deposition or settling function [kg/m2/s]

The entrainment and deposition functions are calculated as

Eb=εscz|z=a=ωsca* (2-44)
Eb=ωsca(245)

where

z = vertical coordinate from the bed [m]
εs = vertical sediment mixing coefficient [m2/s]
c = local sediment concentration [kg/m3]
ωs = sediment fall velocity [m/s]
ca = calculated sediment concentration at an elevation a above the bed [kg/m3]
ca* = equilibrium (capacity) sediment concentration at an elevation a above the bed [kg/m3]

Bed Change Equation

The bed change is calculated as

ρs(1pm')fmorphzbt=qb*jxj+xj(Dsqt*zbxj)+EbDb (2-46)

where

zb = bed elevation with respect to the vertical datum [m]
pm' = bed porosity [-]
fmorph = morphologic acceleration factor [-]
ρs = sediment density [~2650 kg/m3 for quartz sediment]
Ds = empirical bed-slope coefficient (constant) [-]
qbk=hUCtk(1rsk) is the bed load mass transport rate [kg/m/s]

As mentioned above the equilibrium total-load sediment transport model is a single-size sediment transport model and is only available with the explicit time-stepping scheme. For more information on the equilibrium sediment transport model, the reader is referred to Buttolph et al. (2007).

Nonequilibrium Total-Load Transport Model

Total-load Transport Equation

The single-sized sediment transport model described in Sánchez and Wu (2011a) was extended to multiple-sized sediments within CMS by Sánchez and Wu (2011b). In this model, the sediment transport is separated into current- and wave-related transports. The transport due to currents includes the stirring effect of waves; and the wave-related transport includes the transport due to asymmetric oscillatory wave motion as well as steady contributions by Stokes drift, surface roller, and undertow. The current-related bed and suspended transports are combined into a single total-load transport equation, thus reducing the computational costs and simplifying the bed change computation. The two-dimensional horizontal (2DH) transport equation for the current-related total load is

t(hCtkβtk)+(hVjCtk)xj=xj[vsh(rskCtk)xj]+αtωsk(Ct*kCtk)(247)

for j=1,2; k=1,2,......N , where N is the number of sediment size classes and

t = time [s]
h = water depth [m]
xj = Cartesian coordinate in the jth direction [m]
ctk = depth-averaged total-load sediment mass concentration for size class k defined as ctk=qtk/(Uh) in which qtk is the total-load mass transport [kg/m3]
ctk* = depth-averaged total-load sediment mass concentration for size class k and described in the Equilibrium Concentration and Transport Rates section [kg/m3]
βtk = total-load correction factor described in the Total-Load Cor-rection Factor section [-]
rsk = fraction of suspended load in total load for size class k and is described in Fraction of Suspended Sediments section [-]
vs = horizontal sediment mixing coefficient described in the Hori-zontal Sediment Mixing Coefficient section [m2/s]
αt = total-load adaptation coefficient described in the Adaptation Coefficient section [-]
ωsk = sediment fall velocity [m/s]

The above equation may be applied to single-sized sediment transport by using a single sediment size class (i.e. N=1). The bed composition, however, does not vary when using a single sediment size class. The sediment mass concentrations are used rather than volume concentrations in order to avoid precision errors at low concentrations.

Fraction of Suspended Sediments

In order to solve the system of equations for sediment transport implicitly, the fraction of suspended sediments must be determined explicitly. This is done by assuming

rsk=qskqtkqsk*qtk* (2-48)

where qsk and qtk are the actual fraction of suspended- and total-load transport rates and qsk* and qtk* are the equilibrium fraction of suspended- and total-load transport rates.

Adaptation Coefficient

The total-load adaptation coefficient, αt , is an important parameter in the sediment transport model. There are many variations of this parameter in literature (e.g. Lin 1984, Gallappatti and Vreugdenhil 1985, and Armanini and di Silvio 1986). CMS uses a total-load adaptation coefficient αt that is related to the total-load adaptation length Lt and time Tt by


Lt=Uhαtωs=UTt (2-49)

where

ωs = sediment fall velocity corresponding to the transport grain size for single-sized sediment transport or the median grain size for multiple-sized sediment transport [m/s]
U = depth-averaged current velocity [m/s]
h = water depth [m]

The adaptation length (time) is a characteristic distance (time) for sedi-ment to adjust from non-equilibrium to equilibrium transport. Because the total load is a combination of the bed and suspended loads, the associated adaptation length may be calculated as Lt=rsLs+(1rs)Lb or Lt=max(Ls,Lb) , where Ls and Lb are the suspended- and bed-load adaptation lengths. Ls is defined as

Ls=Uhαωs=UTs (2-50)

in which α and Ts are the adaptation coefficient lengths for suspended load for the adaptation coefficient α can be calculated either empirically or based on analytical solutions to the pure vertical convection-diffusion equation of suspended sediment. One example of an empirical formula is that proposed by Lin (1984)


α=3.25+0.55ln(ωsκu*) (2-51)

where u* is the bed shear stress, and κ is the von Karman constant. Armanini and di Silvio (1986) proposed an analytical equation

1α=δh+(1h)exp[1.5(δh)1/6ωsu*] (2-52)

where δ is the thickness of the bottom layer defined by δ=33z0 and z0 is the zero-velocity distance from the bed. Gallapatti (1983) proposed the following equation to determine the suspended load adaptation time

Ts=hu*exp[(1.5720.12ur)ω*3+(326.832ur2.20470.2)ω*2+(0.1385lnur6.4061)ω*+(0.5467ur2.1963)] (2-53)

where u*c is the current related bottom shear velocity, ur=u*/U , and ω*=ωs/u* .

The bed-load adaptation length, Lb , is generally related to the dimension of bed forms such as sand dunes. Large bed forms are generally proportional to the water depth and therefore the bed load adaptation length can be estimated as Lb=abh in which ab is an empirical coefficient on the order of 5-10. Fang (2003) found that Lb of approximately two or three times the grid resolution works well for field applications. Although limited guidance exists on methods to estimate Lb, the determination of Lb is still empirical and in the developmental stage. For a detailed discussion of the adaptation length, the reader is referred to Wu (2007). In general, it is recommended that the adaptation length be calibrated with field data in order to achieve the best and most reliable results.

Total-Load Correction Factor

The correction factor, βtk , accounts for the vertical distribution of the suspended sediment concentration and velocity profiles, as well as the fact that bed load travels a slower velocity than the depth-averaged current velocity (see Figure 2.3). By definition, βtk is the ratio of the depth-averaged total-load and flow velocities.

fig_2_3.png

Figure 2.3. Schematic of sediment and current vertical profiles.

In a combined bed load and suspended load model, the correction factor is given by

βtk=1rsk/βsk+(1rsk)U/ubk(2-54)

where ubk is the bed load velocity and βsk is the suspended load correction factor and is defined as the ratio of the depth-averaged sediment and flow velocities. Since most sediment is transported near the bed, both the total and suspended load correction factors (βtk and βsk ) are usually less than 1 and typically in the range of 0.3 and 0.7, respectively. By assuming logarithmic current velocity and exponential suspended sediment concentration profiles, an explicit expression for the suspended load correction factor βsk may be obtained as (Sánchez and Wu 2011)

βsk=ahuckdzUahckdz=E1(ϕkA)E1(ϕk)+ln(A/Z)eϕkAln(1/Z)eϕkeϕkA[ln(1/Z)1][1eϕk(1A)] (2-55)

where ϕk=ωskh/ϵ,A=a/h,z=za/h, in which ωsk is the sediment fall velocity for size class k, ϵ is the vertical mixing coefficient, a is a reference height for the suspended load, h is the total water depth, za is the apparent roughness length, and E1 is the exponential integral. The equation can be further simplified by assuming that the reference height is proportional to the roughness height (e.g.a=30za ), so that βsk=βsk(z,ϕk). Figure 2.4 shows a comparison of the suspended load correction factor based on the logarithmic velocity with exponential and Rouse suspended sediment concentration profiles.

Figure 2.4. Suspend load correction factors based on the logarithmic velocity profile and (a) exponential and (b) Rouse suspended sediment profile. The Rouse number is r=ωs/(κu*).

The bed load velocity, ubk , is calculated using the van Rijn (1984a) formula with re-calibrated coefficients from Wu et al. (2006)

ubk=1.64(τb'τcrk1)0.5(s1)gdk (2-56)

where s is the specific gravity, g is the gravitational constant, d50 is the median grain size diameter, τb' is the bed shear stress related to the grain roughness and is determined by τb'=(n'/n)3/2τb where n'=d501/6/20 is the Manning’s coefficient corresponding to the grain roughness and τcrk is the critical bed shear stress.

Bed Change Equation

The fractional bed change is calculated as

ρs(1ρm')fmorph(zbt)k=αtωsk(CtkCtk*)+xj(Dsqbkzbxj) (2-57)

where

zb = bed elevation with respect to the vertical datum [m]
pm' = bed porosity [-]
fmorph = morphologic acceleration factor [-]
ρs = sediment density [~2650 kg/m3 for quartz sediment]
Ctk = depth-averaged total-load sediment mass concentration for size class k defined as Ctk=qtk/(Uh) in which qtk is the total-load mass transport [kg/m3]
Ctk* = depth-averaged total-load equilibrium sediment mass concentration for size class k and described in the Equilibrium Concentration and Transport Rates section [kg/m3]
αt = total-load adaptation coefficient described in the Adaptation Coefficient section [-]
ωsk = sediment fall velocity [m/s]
Ds = empirical bed-slope coefficient (constant) [-]
qbk=hUCtk(1rsk) is the bed load mass transport rate [kg/m/s]

The sediment density is required in the previous equation (2-51) since mass concentrations are used. For a detailed derivation of the above equation, the reader is referred to Sanchez and Wu (2011a). The total bed change is calculated as the sum of the bed change for all size classes

zbt=k(zbt)k (2-58)

The purpose of the morphologic acceleration factor fmorph is to speed-up the bed change so that the simulation time tsim represents approximately the change that would occur in tmorph=fmorphtsim . This factor should be used with caution and only for idealized cases or time periods which are periodic (mainly tidal). If time-varying winds or waves are important pro-cesses for driving sediment transport, then it is recommended to use re-duced or representative wind and wave conditions. Since the CMS runs relatively fast, it is generally recommended to not use the morphologic acceleration factor when validating the sediment transport model using hindcast measurements. If good initial and boundary conditions are available should be available and therefore it. The morphologic acceleration factor is useful however when simulating idealized cases or analyzing project alternatives.

Bed material sorting and layering

Bed sorting is the process in which the bed material changes composition (fraction of each grain size class). The bed is descritized into multiple lay-ers each with a uniform bed composition. The fraction of each size class is then calculated and stored in each layer. The sorting of sediments is then calculating using the mixing or active layer concept (Hirano 1971; Karim and Kennedy 1982; and Wu 1991). The active layer is the top layer of the bed which exchanges material directly with the sediment transport.

The temporal variation of the bed-material gradation in the first (mixing or active) layer is calculated as (Wu 2007)

(δ1p1k)t=(zbt)k+p1k*(δ1tzbt) (2-59)

where δ1 is the thickness of the first layer. p1k* is equal to p1k for zb/tδ1/t0 , and equal to the bed material gradation in the second sediment layer for zb/tδ1/t<0 . The bed-material sorting in the se-cond layer is calculated as

(δ2p2k)t=p1k*(δ1tzbt) (2-60)

where δ2 is the thickness of the second layer, and <\math>p_{2k}</math> is the fraction of the kth sediment size in the second layer. It is noted that there is no material exchanged between the sediment layers below the second layer.

The sediment transport, bed change, and bed gradation equations are solved simultaneously (coupled), but are decoupled from the flow calculation at a given time step. To illustrate the bed layering process, Figure 2.5 shows an example of the temporal evolution of 7 bed layers during erosional and depositional regimes.

Figure 2.5. Schematic showing an example bed layer evolution. Colors indicate layer number and not bed composition.

Mixing Layer Thickness

The mixing layer thickness is calculated as

δ1=min[max(δmin,2d50,0.5Δ),δmax]

where Δ is the beform height, and δmin and δmax are user-specified mini-mum and maximum layer thicknesses, respectively. At the beginning of each time step, the mixing layer thickness is calculated. For cell with a hard (non-erodable) bottom, the mixing layer is calculated as

δ1,hb=min(δ1,zbzhb)(2-62)

where zhb is the elevation of the hard bottom. A hard bottom is a nonerodable bed surface such as bed rock or a coastal structure.

Equilibrium Concentrations and Transport Rates

In order to close the system of equations describing the sediment transport, bed change, and bed sorting equations, the fractional equilibri-um depth-averaged total-load concentration Ctk* must be estimated from an empirical formula. The depth-averaged equilibrium concentration is defined as

Ctk*=qtk*Uh(2-63)

where qtk* is the total-load transport for the kth sediment size class esti-mated from an empirical formula.

For convenience, Ctk* is written in general form as

Ctk*=p1kCtk*(2-64)

where p1k is the fraction of the sediment size k in the first (top) bed layer and Ctk* is the potential equilibrium total-load concentration. The potential concentration Ctk* can be interpreted as the equilibrium concentration for uniform sediment of size dk . Equation 2-57 above is essential for the coupling of sediment transport, bed change, and bed sorting equations.

Equilibrium Transport and Concentration Formulas

Lund-CIRP

Camenen and Larson (2005, 2007, 2008) developed general sediment transport formulas for bed and suspended loads under combined waves and currents. These are referred to as the Lund-CIRP transport formulas. The general transport formulas can be used for both symmetric and asymmetric waves, but for simplicity the waves are assumed to be symmetric. The current-related bed- and suspended-load transport with wave stirring is given by

qb*(s1)gd503=fbρs12ΘcΘcw,mexp(4.5ΘcrΘcw)(2-65)
qb*(s1)gd503=fsρsCRUϵωs[1exp(ωshϵ)] (2-66)

where

qb* = potential equilibrium bed load transport [kg/m/s]
qs* = potential suspended load transport [kg/m/s]
d50 = median grain size [m]
s = sediment specific gravity or relative density [-]
g = gravitational constant (9.81 m/s2)
ρs = sediment density (~2650 kg/m3)
ΘC = Shields parameters due to currents [-]
Θcw,m = mean Shields parameters due to waves and currents [-]
Θcw = maximum Shields parameters due to waves and currents [-]
Θcr = critical Shields parameter [-]
ϵ = vertical sediment diffusivity [m2/s]
cR = reference bed concentration [kg/m3]
fb = bed-load scaling factor (default 1.0) [-]
fs = suspended-load scaling factor (default 1.0) [-]

The reference concentration is given by

cR=ρsAcRΘCW,mexp(4.5ΘcrΘcw) (2-67)

where the coefficient AcR is determined by the following relationship

AcR=0.0035exp(0.3d*) (2-68)

The current-related shear stress is calculated as

τc=ρcbU2 (2-69)

where ρ is the water density,cb is the bed friction coefficient, and U is the current velocity magnitude. The drag coefficient is calculated as

cb=[κln(h/z0)1]2(2-70)

where κ is the von Karman constant (0.4), h is the total water depth, and z0 is the roughness length calculated as z0=ks/30 , where ks is the total bed Nikuradse roughness. The total bed roughness is assumed to be a linear sum of the grain-related roughness ks,d, form-drag (ripple) roughness ks,r, and sediment-related roughness ks,s . Bed forms are also separated into current and wave-related bed forms. The current- and wave-related total roughness is then

ks,cw=ksg+ksr,cw+kss,cw(2-71)

where the subscript c|w indicates either the current (c) or wave (w) related component. The grain-related roughness is estimated as ksg=2d50.

The ripple roughness is calculated as (Soulsby 1997)

ksr,cw=7.5Hr,cw2Lr,cw (2-72)

where Hr,cw and Lr,cw are either the current- or wave-related ripple height and length respectively. The current-related ripple height and length are calculated as


Lr,c=1000d50 (2-73)
Hr,c=Lr,c/7(2-74)

The wave-related ripple height and length are calculated using the expressions proposed by van Rijn (1984b, 1989)

Hr,c={0.22Awforψw<102.8×1013(250ψw)5Awfor10ψw<2500for250ψw(2-75)
Lr,w={1.25Awforψw<101.4×106(250ψw)2.5Awfor10ψw<2500for250ψw(2-76)


where AW is the semi-orbital excursion and ψW is the wave mobility parameter. The semi-orbital excursion is defined as

AW=uwT2π<2-77)

in which uw is the peak bottom orbital velocity and T is the wave period. For random waves, uw=urms and T=Tp are used as representative values. The wave mobility parameter,ψw , is defined as


ψW=uW2(s1)gd50(2-78)

For the Lund-CIRP sediment transport equations, a wave bottom shear stress is calculated as

τW=12ρfwuw2(2-79)

where fw is the wave friction factor calculated using the expression of Swart (1976)

fw={exp(5.21r0.196.0)forr>1.570.3forr1.57(2-80)

where r is the relative roughness defined as r=AW/ksg|.

The current- and wave-related sediment roughnesses are estimated as

kss,cw=5d50Θcw(2-81)

Equation 2-73 above must be solved simultaneously with the expressions for the bottom shear stress because the roughness depends on the stress. The exact solution is approximated using explicit polynomial fits in order to avoid time-consuming iterations in calculating the bed shear stress.

The critical Shields parameter is estimated using the formula proposed by Soulsby (1997)

Θcr=0.31+1.2d*+0.055[1exp(0.02d*)] (2-82)

where d* is the dimensionless grain size

d*=d50[(s1)gv2]1/2(2-83)


where v is the kinematic viscosity.

The sediment fall velocity is calculated using the formula by Soulsby (1997)

ωs=vd[(10.362+1.049d*3)1/210.36](2-84)

where d is the grain size. The vertical sediment diffusivity is calculated as

ϵ=h(Deρ)1/3(2-85)

where De is the total effective dissipation given by

De=kb3Db+kc3Dc+kw3Dw(2-86)

in which kb,kc,andkware coefficients, Db is the wave breaking dissipation, and DcandDw are the bottom friction dissipation due to currents and waves, respectively. The dissipation from bottom friction due to current, Dc, and the dissipation from bottom friction due to waves, Dw , are expressed as

Dcw=τcwu*cw(2-87)

where the subscript c|w indicates either the current (c) or wave (w) relat-ed component, and u*cw and τcw are the current- or wave-related bed shear velocity and stress, respectively. The coefficient kb =0.017, and kcandkw are function of the Schmidt number:

kcw=κ6σcw(2-88)

where σcw is either the current or wave-related Schmidt number calculated from the following relationships (Camenen and Larson 2008):

σc={0.7+3.6sin2.5(πωs2u*c)forωsu*c11+3.3sin2.5(πu*c2ωs)forωsu*c>1

(2-89)


σw={0.09+1.4sin2.5(πωs2u*w)forωsu*w11+0.49sin2.5(πu*w2ωs)forωsu*w>1

(2-90)

In the case of coexisting or combined waves and currents, the wave-current Schmidt number is estimated as

σcw=Xv5σc+(1Xv5)σw

(2-91)

where Xv is a weighting factor equal to Xv=U/(U+uw) .

For multiple-sized (nonuniform) sediments, the fractional equilibrium sediment transport rates are calculated as


qbk*(s1)gdk3=fbζkζ1pbkρs12ΘcΘcw,mexp(4.5ΘcrkΘcw)

(2-92)

qbk*(s1)gdk3=fsζkζ1pbkρscRkUϵkωsk[1exp(ωskhϵ)]

(2-93)

where the subscript k indicates variables which are calculated based only on the sediment size class k. ζk is the hiding and exposure coefficient de-scribed in Hiding and Exposure.

Van Rijn

The van Rijn (1984a,b) sediment transport equations for bed load and suspended load are used with the recalibrated coefficients of van Rijn (2007a,b), as given by


qb*=fbρs0.015Uh(UeUcr(s1)gd50)1.5(d50h)1.2(2-94)
qs*=fsρs0.012Ud50(UeUcr(s1)gd50)2.4d*0.6 (2-95)

where

qb* = equilibrium bed-load transport [kg/m/s]
qs* = equilibrium suspended-load transport [kg/m/s]
d50 = median grain size [m]
s = sediment specific gravity or relative density [-]
g = gravitational constant (9.81 m/2)
ρs = sediment density (~2650 kg/m3)
U = depth-averaged current velocity [m/s]
Ucr = critical depth-averaged velocity for incipiant motion [m/s]
Ue = effective depth averaged velocity [m/s]
fs = suspended-load scaling factor (default 1.0) [-]
fb = bed-load scaling factor (default 1.0) [-]

The effective depth-averaged velocity is calculated as Ue=U+γuw with γ =0.4 for random waves and γ =0.8 for regular waves. Here uw is the bottom wave orbital velocity based on linear wave theory and the signifi-cant wave height The critical depth-averaged velocity is estimated as Ucr=βUcrc+(1β)ucrw where β is a blending factor and Ucrc and ucrw are the critical velocities for currents and waves, respectively. As described in van Rijn (2007), the critical velocity for currents is

Ucr={0.19d500.1log10(4hd90),for0.1d500.5mm8.5d500.6log10(4hd90),for0.5d502.0mm(2-96)

where d90 is the sediment grain size in meters of 90th percentile. The critical velocity for waves is based on Komar and Miller (1975),

Ucrw={0.24[(s1)g]0.66d500.33Tp0.33,for0.1d500.5mm0.95[(s1)g]0.57d500.43Tp0.14,for0.5d502.0mm(2-97)

where Tp is the peak wave period. According to van Rijn (2007), the bed load transport formula predicts transport rates by a factor of 2 for velocities higher than 0.6 m/s, but underpredicts transport rates by a factor of 2-3 for velocities close to the initiation of motion.

The near-bed equilibrium sediment concentration for the van Rijn (formula is given by (van Rijn 1985)

ca*=fsρs0.015d50a(τmax'τcrτcr)1.5d*0.3(2-98)

where rmax1 is the maximum skin shear stress, rcr is the critical shear stress, and α is reference height above the bed given by

α=max(0.5Hr,0.01h)(2-99)

where Hr is the bed ripple height.

The van Rijn formula (1984a,b, 2007a,b) were originally proposed for well-sorted sediments. The sediment availability is included by multiplication of transport rateswith the fraction of the sediment size class in the upper bed layer. The hiding and exposure is considered by in a correction factor which multiples to the critical velocity. When applied to multiple-sized sediments, the fractional equilibrium transport rate is calculated as

qbk*=fbρsp1k0.015Uh(Ueζk1/2Ucrk(s1)gdk)1.5(dkh)1.2(2-100)
qsk*=fsρsp1k0.012Uh(Ueζk1/2Ucrk(s1)gdk)2.4(dkh)d*k0.6 (2-101)

where pbk is the fractional bed composition and ζk is the hiding and ex-posure coefficient. The subscript k indicates values which are calculated based on the kth sediment size class.

Soulsby-van Rijn

Soulsby (1997) proposed the following equation for the total load sediment transport rate under action of combined current and waves,

qb*=fbρs0.005Uh(UeUcr(s1)gd50)2.4(d50h)1.2(2-102)
qs*=fbρs0.012Uh(UeUcr(s1)gd50)2.4(d50h)d*0.6(2-103)

where

qt* = equilibrium total-load transport [kg/m/s]
ρs = sediment density (~2650 kg/m3)
U = depth-averaged current velocity [m/s]
Ue = effective velocity defined as =U2+0.018Cdurms2 [m/s]
urms = peak bottom wave orbital velocity based on the root-mean-squared wave height [m/s]
Cd = drag coefficient due to currents alone [-]
Ucr = critical depth-averaged velocity for initiation of motion [m/s]
αsv = empirical bed-slope coefficient, approximately equal to 1.6
ms = bed-slope [-]

The coefficients Asb and Ass are related to the bed and suspended loads, respectively, and given by

Asb=fb0.005h(d50/h)1.2[(s1)gd50]1.2,Ass=fs0.012d50d*0.6[(s1)gd50]1.2(2-104)

Note that the bed and suspended load scaling factors (fb and fs ) are in-cluded in the calculations for the coefficients Asb and Ass.

The current drag coefficient is calculated as

Cd=[κln(h/z0)1]2(2-105)

where z0 is the bed roughness length (equal to 0.006 m).

The Soulsby-van Rijn formula is modified for multiple-sized sediments similarly to the van Rijn formula in the previous section with the equation

qbk*=fbρsp1k0.005Uh(Ueζk1/2Ucrk(s1)gdk)2.4(dkh)1.2(1αsvtanβ)(2-106)
qsk*=fbρsp1k0.012Uh(Ueζk1/2Ucrk(s1)gdk)2.4(dkh)d*k0.6(1esωstanβ)(2-107)

where

Ue=U2+0.018Cdurms2

The subscript k indicates that the value is calculated based only on the size class k and not the median grain size. The availability of sediment fractions is included through p1k, while hiding and exposure of grain sizes is accounted for by modifying the critical velocity. it is noted that the Soulsby-van Rijn (1997) formula are very similar to the van Rijn (1984ab, 207ab) except for the definition of the effective velocity and the recalibration of the bed-load formula coefficients in 2007a. The proposed changes should be verified with measurements or numerical simulations for multiple-sized sediment transport.

Watanabe

The equilibrium total load sediment transport rate is determined by Watanabe (1987) as

qt*=AWatU(τbmaxτcrρg) (2-108)

where

qt* = equilibrium total-load transport [kg/m/s]
τbmax = combined wave-current maximum shear stress [Pa]
ρ = water density (~1025 kg/m3)
g = gravitational constant (9.81 m/s2)
τcr = critical shear stress of incipient motion in [Pa]
AWat = empirical coefficient typically ranging from 0.1 to 2

The critical shear stress is determined from the Shields diagram. The maximum bed shear stress τbmax is calculated as (Soulsby 1997)

τbmax=(τb+τwcosφ)2+(rwsinφ)2(2-109)

where τb is the mean shear stress by waves and current over a wave cycle, τw is the mean wave bed shear stress, and φ is the angle between the waves and current. The wave bed shear stress is calculated as

τw=12ρfwuw2(2-110)

where fw is a wave friction factor and uw is the bottom wave orbital velocity amplitude. The wave friction factor is calculated with the expression by Nielson (1992)

fw=exp(5.5r0.26.3)(2-111)

where r is the relative roughness defined as r=Aw/ksg. Aw is the semi-orbital excursion defined as Aw=uwT/(2π).

Hiding and Exposure

When the bed material is composed of multiple grain sizes, larger grains have a greater probability of being exposed to the flow while smaller particles have a greater probability of being hidden to the flow. Figure 2.6 shows an example of a sediment grain dk being exposed to the flow by an exposure height Δe , and sediment grain dj being hidden by dk.

Figure 2.6. Schematic of the exposure height of bed sediment grains

For the van Rijn (Reference(s)), Soulsby-van Rijn (Reference) and Watanabe (Reference) transport formulas, the hiding and exposure mechanism is considered by correcting the critical shear stress or velocity using a hiding and exposure correction function,ζk . For the Lund-CIRP transport formula, an alternate approach is required due to the way in which the Shields number and grain size are included in the formulation; thus, the hiding and exposure correction function is directly used to multiply the transport rate. Two methods are used to calculate ζk , depending on whether the sediment transport model is run with a single sediment size or with multiple sediment sizes and the methods are described in the following sections.

Single-sized sediment transport

In some applications, the coastal bed material is dominated by a single sediment size with patches of other sediment sizes or materials (e.g. shell hash) that may not contribute significantly to morphology change in the areas of interest; however, they may modify the sediment transport through hiding and exposure. For example, it is possible for the bed material to consist of mostly uniform sand with patches of shell fragments (bimodal distribution) in some regions. The shell material is difficult to model numerically because it is usually poorly sorted and its hydraulic properties are unknown. For such regions, sediment transport models commonly estimate excessive erosion due to the absence of the hiding effect of the coarser shell material (e.g. Cayocca 2001). A better and more physical plausible approach is to use the local bed composition along with a correction function for hiding and exposure to account for the effects of the uniform sand with the patches of coarser shell material. For single-sized sediment transport, the correction function for hiding and exposure is calculated following Parker et al. (1982) and others as,

ζk=(d50dk)(2-112)

where i is an empirical coefficient between 0.5-1.0. The aforementioned sediment transport equations are implemented by using the transport grain size dk rather than the bed material d50 . A single and constant transport size dk is used, while the bed material d50 varies spatially. The spatial distribution of can be obtained from field measurement data and for simplicity is assumed constant during the model simulation time. This is a significant assumption and may not be reasonable for some applications. However, this method provides a simple conceptual mechanism for considering an important process in the proposed single-sized sediment transport model. The approach has been successfully applied to Shinnecock Inlet, NY to simulate morphology change at a coastal inlet (Sánchez and Wu 2011a). A more accurate and complex approach is to simulate the transport and sorting of multiple sediments.

Multiple-sized sediment transport

The hiding and exposure of the each sediment size class is considered by modifying the critical shields parameter Θck for each sediment size class based on Wu et al. 2000

ζk=(PekPhk)m(2-113)

where m is an empirical coefficient that varies for each transport formula, approximately equal to 0.6-1.0. Pek and Phk are the total hiding and exposure probabilities and are calculated as

Phk=j=1Np1jdjlk+djPek=j=1Np1jdkdk+dj (2-114)


where N is the number of grain size classes.

Horizontal Sediment Mixing Coefficient

The horizontal sediment mixing coefficient, vs, represents the combined effects of turbulent diffusion and dispersion due to nonuniform vertical profiles. In CMS, the horizontal sediment mixing coefficient is assumed to be proportional to the turbulent eddy viscosity as

vs=vt/σs(2-115)

where σs is the Schmidt number and vt is…. There are many formulas to estimate the Schmidt number. The Schmidt number is set by default equal to 1.0 in CMS, but may be modified by the user with cards.

Boundary Conditions

At the interface between wet and dry cells, the sediment transport rate is set to zero. The inflow boundary condition requires a given sediment concentration at the boundary. However, for most coastal applications, the actual sediment concentration is not available and the model implements the equilibrium concentration capacity (Dirichlet boundary condition). CMS also requires the size distribution of the inflow transport. For stability reasons, the inflow size distribution is assumed to be equal to the initial size distribution of the bed at the boundary. Inflow sediment transport rates are specified either as a total sediment transport rate, Qsed,tot in kg/s, or as a fractional sediment transport rate,Qsed,k in kg/s, for each sediment size class. If only the total sediment transport rate is specified, then the fractional sediment transport rate is calculated as Qsed,k=p1kQsed.tot where p1k is the bed material fraction of the first layer at the boundary. The fractional sediment transport qtk,B at boundary cell B in kg/m/s is then calculated along the cell string according to

qtk,B=qtk,B*ΔlfQtk,N*Qsed,k(2-116)

where qtk,B* is the potential sediment transport rate at boundary cell B and Δlf is the inflow cell face width. If the flow is directed outward of the do-main, a zero-gradient boundary condition is used for sediment concentration (Neumann boundary condition).

Salinity Transport

Transport Equation

CMS calculates the salinity transport based on the following 2DH salinity conservation equation

(hCsal)t+(hVjCsal)xj=xf(vsalhCsalxj)+Ssal(2-117)

where

Csal = depth-averaged salinity [ppt or parts per thousand]
h = water depth [m]
Vj = total flux velocity [m/s]
vsal = horizontal mixing coefficient vsal=vt/σs[m2/s]
Ssal = source/sink term due to precipitation, evaporation, and structures (e.g. culverts) [ppt m/s]

Equation (2-106) represents the horizontal fluxes of salt in water bodies and is balanced by exchanges of salt via diffusive fluxes. Major processes that contribute to the salinity are: freshwater inflows from rivers, vertical fluxes of freshwater by precipitation and evaporation at the water surface, and groundwater fluxes. These processes can be specified as the surface and bottom boundary conditions in the equation.

Boundary Conditions

The boundary condition for the salinity transport equation is dependent on the flow direction. At cell faces between wet and dry cells, a zero-flux boundary condition is applied. The cell face velocity Uf is zero between wet and dry cells and this eliminates advection transport. In addition, the salinity horizontal mixing coefficient at the cell face, v¯sal,f , is set to zero and this eliminates the diffusive flux. If the flow is directed inward of the modeling domain, then a user-specified salinity concentration is pre-scribed (Direchlet boundary condition). If the flow is directed outward of the modeling domain, then a zero-gradient boundary condition is applied (Neumann boundary condition). Overall comments for this section: 1) A figure would be very HELPFUL showing inward and outward directed flow with the modeling domain and 2) Direchlet and Neumann BCs should be explained with reference citations.

Surface Roller

As a wave transitions from non-breaking to fully-breaking, some of the energy is converted into momentum that goes into the aerated region of the water column known. This phenomenon is known as the surface roller. Under the assumption that the surface roller moves in the mean wave direction, the evolution and dissipation of the surface roller energy is calculated by a steady-steate energy balance equation (Stive and De Vriend 1994, Ruessink et al. 2001)

(2Esrcwj)xj=Dsr+feDbr(2-118)

where

Esr = surface roller energy density [N/m]
c = roller propagation speed [m/s]
wj=(cosθm,sinθm) is the wave unit vector [-]
θm = mean wave direction [deg]
Dsr = roller energy dissipation [N/m/s]
Dbr = wave breaking dissipation (from the wave model) [N/m/s]

The surface roller speed is calculated using the long-wave approximation c=gh , where h is the local water depth. The surface roller dissipation is approximated as

Dsr=g2EsrβDc(2-119)

where βD is the surface roller dissipation coefficient approximately equal to 0.05-0.1. The surface roller contribution to the wave stresses, Rij , is given by

Rij=2Esrwiwj(2-120)

One effect of the surface roller is to push the peak alongshore current ve-locity closer to shore. Another side effect of the surface roller is to improve model stability (Sánchez et al. 2011a). The influence of the surface roller on the mean water surface elevation is relatively minor (Sánchez et al. 2011a).