Sediment Transport 1

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: $q_{tk*}=q_{sk*}+q_{bk*}$.

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.

${\frac {\rho _{s}(1-\rho _{m}^{'})}{f_{morph}}}{\frac {\partial z_{b}}{\partial t}}={\frac {\partial q_{t*j}}{\partial x_{j}}}+{\frac {\partial }{\partial x_{j}}}\left(D_{s}\mid q_{t*}\mid {\frac {\partial z_{b}}{\partial x_{j}}}\right)$ (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]

$x_{j}$ =Cartesian coordinate in the j^th direction [m]

$q_{t*}$ = equilibrium total-load transport rate [kg/m/s]

$z_{b}$ = bed elevation with respect to the vertical datum [m]

$p_{m}^{'}$ = bed porosity [-]

$f_{morph}$ = morphologic acceleration factor [-]

$\rho _{s}$ = sediment density [~2650 kg/m³ for quartz sediment]

$D_{s}$ = 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

${\frac {\partial (hC)}{\partial t}}+{\frac {\partial (hV_{j}C)}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left(v_{s}h{\frac {\partial C}{\partial x_{j}}}\right)+E_{b}-D_{b}$ (2-43)

where

t = time[s]

h = water depth [m]

$x_{j}$ = Cartesian coordinate in the

$E_{b}$ entrainment or pick-up function [kg/m²/s]

$D_{b}$ deposition or settling function [kg/m²/s]

The entrainment and deposition functions are calculated as

$E_{b}=-\varepsilon _{s}{\frac {\partial c}{\partial z}}\left|_{z=a}\right.=\omega _{s}c_{a*}$ (2-44)

$E_{b}=\omega _{s}c_{a}(2-45)$

where

z = vertical coordinate from the bed [m]

$\varepsilon _{s}$ = vertical sediment mixing coefficient [m²/s]

c = local sediment concentration [kg/m³]

$\omega _{s}$ = sediment fall velocity [m/s]

$c_{a}$ = calculated sediment concentration at an elevation a above the bed [kg/m³]

$c_{a*}$ = equilibrium (capacity) sediment concentration at an elevation a above the bed [kg/m³]

Bed Change Equation

The bed change is calculated as

${\frac {\rho _{s}(1-p_{m}^{'})}{f_{morph}}}{\frac {\partial z_{b}}{\partial t}}={\frac {\partial q_{b*j}}{\partial x_{j}}}+{\frac {\partial }{\partial x_{j}}}\left(D_{s}\mid q_{t}*\mid {\frac {\partial z_{b}}{\partial x_{j}}}\right)+E_{b}-D_{b}$ (2-46)

where

$z_{b}$ = bed elevation with respect to the vertical datum [m]

$p_{m}^{'}$ = bed porosity [-]

$f_{morph}$ = morphologic acceleration factor [-]

$\rho _{s}$ = sediment density [~2650 kg/m³ for quartz sediment]

$D_{s}$ = empirical bed-slope coefficient (constant) [-]

$q_{bk}=hUC_{tk}(1-r_{sk})$ 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

${\frac {\partial }{\partial t}}\left({\frac {hC_{tk}}{\beta _{tk}}}\right)+{\frac {\partial (hV_{j}C_{tk})}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left[v_{s}h{\frac {\partial (r_{sk}C_{tk})}{\partial x_{j}}}\right]+\alpha _{t}\omega _{sk}(C_{t*k}-C_{tk})(2-47)$

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]

$x_{j}$ = Cartesian coordinate in the j^th direction [m]

$c{tk}$ = depth-averaged total-load sediment mass concentration for size class k defined as $c_{tk}=q_{tk}/(Uh)$ in which $q_{tk}$ is the total-load mass transport [kg/m³]

$c_{tk*}$ = depth-averaged total-load sediment mass concentration for size class k and described in the Equilibrium Concentration and Transport Rates section [kg/m³]

$\beta _{tk}$ = total-load correction factor described in the Total-Load Cor-rection Factor section [-]

$r_{sk}$ = fraction of suspended load in total load for size class k and is described in Fraction of Suspended Sediments section [-]

$v_{s}$ = horizontal sediment mixing coefficient described in the Hori-zontal Sediment Mixing Coefficient section [m²/s]

$\alpha _{t}$ = total-load adaptation coefficient described in the Adaptation Coefficient section [-]

$\omega _{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

$r_{sk}={\frac {q_{sk}}{q_{tk}}}\Box {\frac {q_{sk*}}{q_{tk*}}}$ (2-48)

where $q_{sk}$ and $q_{tk}$ are the actual fraction of suspended- and total-load transport rates and $q_{sk*}$ and $q_{tk*}$ are the equilibrium fraction of suspended- and total-load transport rates.

Adaptation Coefficient

The total-load adaptation coefficient, $\alpha _{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 $\alpha _{t}$ that is related to the total-load adaptation length $L_{t}$ and time $T_{t}$ by

$L_{t}={\frac {Uh}{\alpha _{t}\omega _{s}}}=UT_{t}$ (2-49)

where

$\omega _{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 $L_{t}=r_{s}L_{s}+(1-r_{s})L_{b}$ or $L_{t}=max(L_{s},L_{b})$ , where L_s and L_b are the suspended- and bed-load adaptation lengths. L_s is defined as

$L_{s}={\frac {Uh}{\alpha \omega _{s}}}=UT_{s}$ (2-50)

in which $\alpha$ and $T_{s}$ are the adaptation coefficient lengths for suspended load for the adaptation coefficient $\alpha$ 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)

$\alpha =3.25+0.55ln\left({\frac {\omega _{s}}{\kappa u_{*}}}\right)$ (2-51)

where $u_{*}$ is the bed shear stress, and $\kappa$ is the von Karman constant. Armanini and di Silvio (1986) proposed an analytical equation

${\frac {1}{\alpha }}={\frac {\delta }{h}}+\left(1-{\frac {\partial }{h}}\right)exp\left[-1.5\left({\frac {\delta }{h}}\right)^{-1/6}{\frac {\omega _{s}}{u_{*}}}\right]$ (2-52)

where $\delta$ is the thickness of the bottom layer defined by $\delta =33z_{0}$ and $z_{0}$ is the zero-velocity distance from the bed. Gallapatti (1983) proposed the following equation to determine the suspended load adaptation time

$T_{s}={\frac {h}{u_{*}}}exp\left[{\begin{aligned}&(1.57-20.12u_{r})\omega _{*}^{3}+(326.832u_{r}^{2.2047}-0.2)\omega _{*}^{2}\\&+(0.1385lnu_{r}-6.4061)\omega _{*}+(0.5467u_{r}-2.1963)\end{aligned}}\right]$ (2-53)

where $u*_{c}$ is the current related bottom shear velocity, $u_{r}=u*/U$ , and $\omega _{*}=\omega _{s}/u_{*}$ .

The bed-load adaptation length, $L_{b}$ , 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 $L_{b}=a_{b}h$ in which $a_{b}$ is an empirical coefficient on the order of 5-10. Fang (2003) found that $L_{b}$ of approximately two or three times the grid resolution works well for field applications. Although limited guidance exists on methods to estimate $L_{b}$, the determination of $L_{b}$ 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, $\beta _{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, $\beta _{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

$\beta _{tk}={\frac {1}{r_{sk}/\beta _{sk}+(1-r_{sk})U/u_{bk}}}$(2-54)

where $u_{bk}$ is the bed load velocity and $\beta _{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 ($\beta _{tk}$ and $\beta _{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 $\beta _{sk}$ may be obtained as (Sánchez and Wu 2011)

$\beta _{sk}={\frac {\int _{a}^{h}uc_{k}dz}{U\int _{a}^{h}c_{k}dz}}={\frac {E_{1}(\phi _{k}A)-E_{1}(\phi _{k})+ln(A/Z)e^{-\phi _{k}A}-ln(1/Z)e^{-\phi _{k}}}{e^{-\phi _{k}A}[ln(1/Z)-1][1-e^{-\phi _{k}(1-A)}]}}$ (2-55)

where $\phi _{k}=\omega _{sk}h/\epsilon ,A=a/h,z=z_{a}/h,$ in which $\omega _{sk}$ is the sediment fall velocity for size class k, $\epsilon$ is the vertical mixing coefficient, a is a reference height for the suspended load, h is the total water depth, $z_{a}$ is the apparent roughness length, and $E_{1}$ 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=30z_{a}$ ), so that $\beta _{sk}=\beta _{sk}(z,\phi _{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=\omega _{s}/(\kappa u_{*})$.

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

$u_{bk}=1.64\left({\frac {\tau _{b}^{'}}{\tau _{crk}}}-1\right)^{0.5}{\sqrt {(s-1)gd_{k}}}$ (2-56)

where s is the specific gravity, g is the gravitational constant, $d_{50}$ is the median grain size diameter, $\tau _{b}^{'}$ is the bed shear stress related to the grain roughness and is determined by $\tau _{b}^{'}=(n^{'}/n)^{3/2}\tau _{b}$ where $n^{'}=d_{50}^{1/6}/20$ is the Manning’s coefficient corresponding to the grain roughness and $\tau _{crk}$ is the critical bed shear stress.

Bed Change Equation

The fractional bed change is calculated as

${\frac {\rho _{s}(1-\rho _{m}^{'})}{f_{morph}}}\left({\frac {\partial z_{b}}{\partial t}}\right)_{k}=\alpha _{t}\omega _{sk}(C_{tk}-C_{tk*})+{\frac {\partial }{\partial x_{j}}}\left(D_{s}\mid q_{bk}\mid {\frac {\partial z_{b}}{\partial x_{j}}}\right)$ (2-57)

where

$z_{b}$ = bed elevation with respect to the vertical datum [m]

$p_{m}^{'}$ = bed porosity [-]

$f_{morph}$ = morphologic acceleration factor [-]

$\rho _{s}$ = sediment density [~2650 kg/m³ for quartz sediment]

$C_{tk}$ = depth-averaged total-load sediment mass concentration for size class k defined as $C_{tk}=q_{tk}/(Uh)$ in which $q_{tk}$ is the total-load mass transport [kg/m³]

$C_{tk*}$ = depth-averaged total-load equilibrium sediment mass concentration for size class k and described in the Equilibrium Concentration and Transport Rates section [kg/m³]

$\alpha _{t}$ = total-load adaptation coefficient described in the Adaptation Coefficient section [-]

$\omega _{sk}$ = sediment fall velocity [m/s]

$D_{s}$ = empirical bed-slope coefficient (constant) [-]

$q_{bk}=hUC_{tk}(1-r_{sk})$ 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

${\frac {\partial z_{b}}{\partial t}}=\sum _{k}\left({\frac {\partial z_{b}}{\partial t}}\right)_{k}$ (2-58)

The purpose of the morphologic acceleration factor $f_{morph}$ is to speed-up the bed change so that the simulation time $t_{sim}$ represents approximately the change that would occur in $t{morph}=f_{morph}t_{sim}$ . 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)

${\frac {\partial (\delta _{1}p_{1k})}{\partial t}}=\left({\frac {\partial z_{b}}{\partial t}}\right)_{k}+p_{1k}^{*}\left({\frac {\partial \delta _{1}}{\partial t}}-{\frac {\partial z_{b}}{\partial t}}\right)$ (2-59)

where $\delta _{1}$ is the thickness of the first layer. $p_{1k}^{*}$ is equal to $p_{1k}$ for $\partial z_{b}/\partial t-\partial \delta _{1}/\partial t\geq 0$ , and equal to the bed material gradation in the second sediment layer for $\partial z_{b}/\partial t-\partial \delta _{1}/\partial t<0$ . The bed-material sorting in the se-cond layer is calculated as

${\frac {\partial (\delta _{2}p_{2k})}{\partial t}}=-p_{1k}^{*}\left({\frac {\partial \delta _{1}}{\partial t}}-{\frac {\partial z_{b}}{\partial t}}\right)$ (2-60)

where $\delta _{2}$ is the thickness of the second layer, and <\math>p_{2k}</math> is the fraction of the $k^{th}$ 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

$\delta _{1}=min[max(\delta _{min},2d_{50},0.5\Delta ),\delta _{max}]$

where $\Delta$ is the beform height, and $\delta _{min}$ and $\delta _{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

$\delta _{1,hb}=min(\delta _{1},z_{b}-z_{hb})$(2-62)

where $z_{hb}$ 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 $C_{tk*}$ must be estimated from an empirical formula. The depth-averaged equilibrium concentration is defined as

$C_{tk*}={\frac {q_{}tk*}{Uh}}$(2-63)

where $q_{tk*}$ is the total-load transport for the $k^{th}$ sediment size class esti-mated from an empirical formula.

For convenience, $C_{tk*}$ is written in general form as

$C_{tk*}=p_{1k}C_{tk}^{*}$(2-64)

where $p_{1k}$ is the fraction of the sediment size k in the first (top) bed layer and $C_{tk}^{*}$ is the potential equilibrium total-load concentration. The potential concentration $C_{tk}^{*}$ can be interpreted as the equilibrium concentration for uniform sediment of size $d_{k}$ . 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

${\frac {q_{b}^{*}}{\sqrt {(s-1)gd_{50}^{3}}}}=f_{b}\rho _{s}12{\sqrt {\Theta _{c}}}\Theta _{cw,m}exp\left(-4.5{\frac {\Theta _{cr}}{\Theta _{cw}}}\right)$(2-65)

${\frac {q_{b}^{*}}{\sqrt {(s-1)gd_{50}^{3}}}}=f_{s}\rho _{s}C_{R}U{\frac {\epsilon }{\omega _{s}}}\left[1-exp\left(-{\frac {\omega _{s}h}{\epsilon }}\right)\right]$ (2-66)

where

$q_{b}^{*}$ = potential equilibrium bed load transport [kg/m/s]

$q_{s}^{*}$ = potential suspended load transport [kg/m/s]

$d_{50}$ = median grain size [m]

$s$ = sediment specific gravity or relative density [-]

$g$ = gravitational constant (9.81 m/s²)

$\rho _{s}$ = sediment density (~2650 kg/m³)

$\Theta _{C}$ = Shields parameters due to currents [-]

$\Theta _{cw,m}$ = mean Shields parameters due to waves and currents [-]

$\Theta _{c}w$ = maximum Shields parameters due to waves and currents [-]

$\Theta _{c}r$ = critical Shields parameter [-]

$\epsilon$ = vertical sediment diffusivity [m²/s]

$c_{R}$ = reference bed concentration [kg/m³]

$f_{b}$ = bed-load scaling factor (default 1.0) [-]

$f_{s}$ = suspended-load scaling factor (default 1.0) [-]

The reference concentration is given by

$c_{R}=\rho _{s}A_{cR}\Theta _{CW,m}exp\left(-4.5{\frac {\Theta _{cr}}{\Theta _{cw}}}\right)$ (2-67)

where the coefficient $A_{cR}$ is determined by the following relationship

$A_{cR}=0.0035exp(-0.3d_{*})$ (2-68)

The current-related shear stress is calculated as

$\tau _{c}=\rho c_{b}U^{2}$ (2-69)

where $\rho$ is the water density,$c_{b}$ is the bed friction coefficient, and $U$ is the current velocity magnitude. The drag coefficient is calculated as

$c_{b}=\left[{\frac {\kappa }{ln(h/z_{0})-1}}\right]^{2}$(2-70)

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

$k_{s,c\mid w}=k_{sg}+k_{sr,c\mid w}+k_{ss,c\mid w}$(2-71)

where the subscript c|w indicates either the current (c) or wave (w) related component. The grain-related roughness is estimated as $k_{sg}=2d_{50}$.

The ripple roughness is calculated as (Soulsby 1997)

$k_{sr,c\mid w}=7.5{\frac {H_{r,c\mid w}^{2}}{L_{r,c\mid w}}}$ (2-72)

where $H_{r,c\mid w}$ and $L_{r,c\mid w}$ are either the current- or wave-related ripple height and length respectively. The current-related ripple height and length are calculated as

$L_{r,c}=1000d_{50}$ (2-73)

$H_{r,c}=L_{r,c}/7$(2-74)

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

$H_{r,c}=\left\{{\begin{aligned}&0.22A_{w}&{\text{for}}\psi _{w}<10\\&2.8\times 10^{-13}(250-\psi _{w})^{5}A_{w}&{\text{for}}10\leq \psi _{w}<250\\&0&{\text{for}}250\leq \psi _{w}\end{aligned}}\right.$(2-75)

$L_{r,w}=\left\{{\begin{aligned}&1.25A_{w}&{\text{for}}\psi _{w}<10\\&1.4\times 10^{-6}(250-\psi _{w})^{2.5}A_{w}&{\text{for}}10\leq \psi _{w}<250\\&0\quad &{\text{for}}250\leq \psi _{w}\end{aligned}}\right.$(2-76)

where $A_{W}$ is the semi-orbital excursion and $\psi _{W}$ is the wave mobility parameter. The semi-orbital excursion is defined as

$A_{W}={\frac {u_{w}T}{2\pi }}$<2-77)

in which $u_{w}$ is the peak bottom orbital velocity and $T$ is the wave period. For random waves, $u_{w}={\sqrt {u}}_{rms}$ and $T=T_{p}$ are used as representative values. The wave mobility parameter,$\psi _{w}$ , is defined as

$\psi _{W}={\frac {u_{W}^{2}}{(s-1)gd_{50}}}$(2-78)

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

$\tau _{W}={\frac {1}{2}}\rho f_{w}u_{w}^{2}$(2-79)

where $f_{w}$ is the wave friction factor calculated using the expression of Swart (1976)

$f_{w}=\left\{{\begin{aligned}&exp(5.21r^{-0.19}-6.0){\text{for}}r>1.57\\&0.3\quad \quad \quad \quad \quad \quad \quad \quad {\text{for}}r\leq 1.57\end{aligned}}\right.$(2-80)

where $r$ is the relative roughness defined as $r=A_{W}/k_{sg|}$.

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

$k_{ss,c\mid w}=5d_{50}\Theta _{c\mid w}$(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)

$\Theta _{cr}={\frac {0.3}{1+1.2d_{*}}}+0.055[1-exp(-0.02d_{*})]$ (2-82)

where $d_{*}$ is the dimensionless grain size

$d_{*}=d_{50}\left[{\frac {(s-1)g}{v_{2}}}\right]^{1/2}$(2-83)

where $v$ is the kinematic viscosity.

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

$\omega _{s}={\frac {v}{d}}\left[(10.36^{2}+1.049d_{*}^{3})^{1/2}-10.36\right]$(2-84)

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

$\epsilon =h\left({\frac {D_{e}}{\rho }}\right)^{1/3}$(2-85)

where $D_{e}$ is the total effective dissipation given by

$D_{e}=k_{b}^{3}D_{b}+k_{c}^{3}D_{c}+k_{w}^{3}D_{w}$(2-86)

in which $k_{b},k_{c},{\text{and}}k_{w}$are coefficients, $D_{b}$ is the wave breaking dissipation, and $D_{c}{\text{and}}D_{w}$ are the bottom friction dissipation due to currents and waves, respectively. The dissipation from bottom friction due to current, $D_{c}$, and the dissipation from bottom friction due to waves, $D_{w}$ , are expressed as

$D_{c\mid w}=\tau _{c\mid w}u*_{c\mid w}$(2-87)

where the subscript c|w indicates either the current (c) or wave (w) relat-ed component, and $u*_{c\mid w}$ and $\tau _{c\mid w}$ are the current- or wave-related bed shear velocity and stress, respectively. The coefficient $k_{b}$ =0.017, and $k_{c}{\text{and}}k_{w}$ are function of the Schmidt number:

$k_{c\mid w}={\frac {\kappa }{6}}\sigma _{c\mid w}$(2-88)

where $\sigma _{c\mid w}$ is either the current or wave-related Schmidt number calculated from the following relationships (Camenen and Larson 2008):

$\sigma _{c}=\left\{{\begin{aligned}&0.7+3.6sin^{2.5}\left({\frac {\pi \omega _{s}}{2u*_{c}}}\right)&{\text{for}}{\frac {\omega _{s}}{u*_{c}}}\leq 1\\&1+3.3sin^{2.5}\left({\frac {\pi u*_{c}}{2\omega _{s}}}\right)&{\text{for}}{\frac {\omega _{s}}{u*_{c}}}>1\end{aligned}}\right.$

(2-89)

$\sigma _{w}=\left\{{\begin{aligned}&0.09+1.4sin^{2.5}\left({\frac {\pi \omega _{s}}{2u*_{w}}}\right)&{\text{for}}{\frac {\omega _{s}}{u*_{w}}}\leq 1\\&1+0.49sin^{2.5}\left({\frac {\pi u*_{w}}{2\omega _{s}}}\right)&{\text{for}}{\frac {\omega _{s}}{u*_{w}}}>1\end{aligned}}\right.$

(2-90)

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

$\sigma _{cw}=X_{v}^{5}\sigma _{c}+(1-X_{v}^{5})\sigma _{w}$

(2-91)

where $X_{v}$ is a weighting factor equal to $X_{v}=U/(U+u_{w})$ .

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

${\frac {q_{bk*}}{\sqrt {(s-1)gd_{k}^{3}}}}=f_{b}\zeta _{k}^{\zeta -1}p_{bk}\rho _{s}12{\sqrt {\Theta _{c}}}\Theta _{cw,m}exp\left(-4.5{\frac {\Theta _{crk}}{\Theta _{cw}}}\right)$

(2-92)

${\frac {q_{bk*}}{\sqrt {(s-1)gd_{k}^{3}}}}=f_{s}\zeta _{k}^{\zeta -1}p{bk}\rho _{s}c_{Rk}U{\frac {\epsilon _{k}}{\omega _{sk}}}\left[1-exp\left(-{\frac {\omega _{sk}h}{\epsilon }}\right)\right]$

(2-93)

where the subscript k indicates variables which are calculated based only on the sediment size class k. $\zeta _{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

$q_{b}^{*}=f_{b}\rho _{s}0.015Uh\left({\frac {U_{e}-U_{cr}}{\sqrt {(s-1)gd_{50}}}}\right)^{1.5}\left({\frac {d_{50}}{h}}\right)^{1.2}$(2-94)

$q_{s}^{*}=f_{s}\rho _{s}0.012Ud_{50}\left({\frac {U_{e}-U_{cr}}{\sqrt {(s-1)gd_{50}}}}\right)^{2.4}d_{*}^{-0.6}$ (2-95)

where

$q_{b}^{*}$ = equilibrium bed-load transport [kg/m/s]

$q_{s}^{*}$ = equilibrium suspended-load transport [kg/m/s]

$d_{50}$ = median grain size [m]

s = sediment specific gravity or relative density [-]

g = gravitational constant (9.81 m/²)

$\rho _{s}$ = sediment density (~2650 kg/m³)

U = depth-averaged current velocity [m/s]

$U_{cr}$ = critical depth-averaged velocity for incipiant motion [m/s]

$U_{e}$ = effective depth averaged velocity [m/s]

$f_{s}$ = suspended-load scaling factor (default 1.0) [-]

$f_{b}$ = bed-load scaling factor (default 1.0) [-]

The effective depth-averaged velocity is calculated as $U_{e}=U+\gamma u_{w}$ with $\gamma$ =0.4 for random waves and $\gamma$ =0.8 for regular waves. Here $u_{w}$ 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 $U_{cr}=\beta U_{crc}+(1-\beta )u_{crw}$ where $\beta$ is a blending factor and $U_{crc}$ and $u_{crw}$ are the critical velocities for currents and waves, respectively. As described in van Rijn (2007), the critical velocity for currents is

$U_{cr}=\left\{{\begin{aligned}&0.19d_{50}^{0.1}log_{10}\left({\frac {4h}{d_{90}}}\right),&for0.1\leq d_{50}\leq 0.5mm\\&8.5d_{50}^{0.6}log_{10}\left({\frac {4h}{d_{90}}}\right),&for0.5\leq d_{50}\leq 2.0mm\end{aligned}}\right.$(2-96)

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

$U_{crw}=\left\{{\begin{aligned}&0.24[(s-1)g]^{0.66}d_{50}^{0.33}T_{p}^{0.33},&for0.1\leq d_{50}\leq 0.5mm\\&0.95[(s-1)g]^{0.57}d_{50}^{0.43}T_{p}^{0.14},&for0.5\leq d_{50}\leq 2.0mm\end{aligned}}\right.$(2-97)

where $T_{p}$ 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)

$c_{a*}=f_{s}\rho _{s}0.015{\frac {d_{50}}{a}}\left({\frac {\tau _{max}^{'}-\tau _{cr}}{\tau _{cr}}}\right)^{1.5}d_{*}^{-0.3}$(2-98)

where $r_{max}^{1}$ is the maximum skin shear stress, $r_{cr}$ is the critical shear stress, and $\alpha$ is reference height above the bed given by

$\alpha =max(0.5H_{r},0.01h)$(2-99)

where $H_{r}$ 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

$q_{bk*}=f_{b}\rho _{s}p_{1k}0.015Uh\left({\frac {U_{e}-\zeta _{k}^{1/2}U_{crk}}{\sqrt {(s-1)gd_{k}}}}\right)^{1.5}\left({\frac {d_{k}}{h}}\right)^{1.2}$(2-100)

$q_{sk*}=f_{s}\rho _{s}p_{1k}0.012Uh\left({\frac {U_{e}-\zeta _{k}^{1/2}U_{crk}}{\sqrt {(s-1)gd_{k}}}}\right)^{2.4}\left({\frac {d_{k}}{h}}\right)d_{*k}^{-0.6}$ (2-101)

where $p_{bk}$ is the fractional bed composition and $\zeta _{k}$ is the hiding and ex-posure coefficient. The subscript k indicates values which are calculated based on the k^th 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,

$q_{b*}=f_{b}\rho _{s}0.005Uh\left({\frac {U_{e}-U_{cr}}{\sqrt {(s-1)gd_{50}}}}\right)^{2.4}\left({\frac {d_{50}}{h}}\right)^{1.2}$(2-102)

$q_{s*}=f_{b}\rho _{s}0.012Uh\left({\frac {U_{e}-U_{cr}}{\sqrt {(s-1)gd_{50}}}}\right)^{2.4}\left({\frac {d_{50}}{h}}\right)d_{*}^{-0.6}$(2-103)

where

$q_{t*}$ = equilibrium total-load transport [kg/m/s]

$\rho _{s}$ = sediment density (~2650 kg/m³)

U = depth-averaged current velocity [m/s]

$U_{e}$ = effective velocity defined as $={\sqrt {U^{2}+{\frac {0.018}{C_{d}}}u_{rms}^{2}}}$ [m/s]

$u{rms}$ = peak bottom wave orbital velocity based on the root-mean-squared wave height [m/s]

$C_{d}$ = drag coefficient due to currents alone [-]

$U_{cr}$ = critical depth-averaged velocity for initiation of motion [m/s]

$\alpha _{sv}$ = empirical bed-slope coefficient, approximately equal to 1.6

$m_{s}$ = bed-slope [-]

The coefficients A_sb and A_ss are related to the bed and suspended loads, respectively, and given by

$A_{sb}=f_{b}{\frac {0.005h(d_{50}/h)^{1.2}}{[(s-1)gd_{50}]^{1.2}}},A_{ss}=f_{s}{\frac {0.012d_{50}d_{*}^{-0.6}}{[(s-1)gd_{50}]^{1.2}}}$(2-104)

Note that the bed and suspended load scaling factors ($f_{b}$ and $f_{s}$ ) are in-cluded in the calculations for the coefficients $A_{sb}$ and $A_{ss}$.

The current drag coefficient is calculated as

$C_{d}=\left[{\frac {\kappa }{ln(h/z_{0})-1}}\right]^{2}$(2-105)

where $z_{0}$ 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

$q_{bk*}=f_{b}\rho _{s}p_{1k}0.005Uh\left({\frac {U_{e}-\zeta _{k}^{1/2}U_{crk}}{\sqrt {(s-1)gd_{k}}}}\right)^{2.4}\left({\frac {d_{k}}{h}}\right)^{1.2}(1-\alpha _{sv}tan\beta )$(2-106)

$q_{sk*}=f_{b}\rho _{s}p_{1k}0.012Uh\left({\frac {U_{e}-\zeta _{k}^{1/2}U_{crk}}{\sqrt {(s-1)gd_{k}}}}\right)^{2.4}\left({\frac {d_{k}}{h}}\right)d_{*k}^{-0.6}\left(1-{\frac {e_{s}}{\omega _{s}}}tan\beta \right)$(2-107)

where

$U_{e}={\sqrt {U^{2}+{\frac {0.018}{C_{d}}}u_{rms}^{2}}}$

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 $p_{1k}$, 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

$q_{t*}=A_{Wat}U\left({\frac {\tau _{bmax}-\tau _{cr}}{\rho g}}\right)$ (2-108)

where

$q_{t*}$ = equilibrium total-load transport [kg/m/s]

$\tau _{bmax}$ = combined wave-current maximum shear stress [Pa]

$\rho$ = water density (~1025 kg/m³)

g = gravitational constant (9.81 m/s²)

$\tau _{cr}$ = critical shear stress of incipient motion in [Pa]

$A_{Wat}$ = empirical coefficient typically ranging from 0.1 to 2

The critical shear stress is determined from the Shields diagram. The maximum bed shear stress $\tau _{bmax}$ is calculated as (Soulsby 1997)

$\tau _{bmax}={\sqrt {(\tau _{b}+\tau _{w}cos\varphi )^{2}+(r_{w}sin\varphi )^{2}}}$(2-109)

where $\tau _{b}$ is the mean shear stress by waves and current over a wave cycle, $\tau _{w}$ is the mean wave bed shear stress, and $\varphi$ is the angle between the waves and current. The wave bed shear stress is calculated as

$\tau _{w}={\frac {1}{2}}\rho f_{w}u_{w}^{2}$(2-110)

where $f_{w}$ is a wave friction factor and $u_{w}$ is the bottom wave orbital velocity amplitude. The wave friction factor is calculated with the expression by Nielson (1992)

$f_{w}=exp(5.5r^{-0.2}-6.3)$(2-111)

where r is the relative roughness defined as $r=A_{w}/k_{s}g$. $A_{w}$ is the semi-orbital excursion defined as $A_{w}=u_{w}T/(2\pi )$.

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 $d_{k}$ being exposed to the flow by an exposure height $\Delta _{e}$ , and sediment grain $d_{j}$ being hidden by $d_{k}$.

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,$\zeta _{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 $\zeta _{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,

$\zeta _{k}=\left({\frac {d_{50}}{d_{k}}}\right)$(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 $d_{k}$ rather than the bed material $d_{50}$ . A single and constant transport size $d_{k}$ is used, while the bed material $d_{50}$ 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 $\Theta _{ck}$ for each sediment size class based on Wu et al. 2000

$\zeta _{k}=\left({\frac {P_{ek}}{P_{hk}}}\right)^{-m}$(2-113)

where m is an empirical coefficient that varies for each transport formula, approximately equal to 0.6-1.0. $P_{ek}$ and $P_{hk}$ are the total hiding and exposure probabilities and are calculated as

$P_{hk}=\sum _{j=1}^{N}p_{1j}{\frac {d_{j}}{l_{k}+d_{j}}}\quad P_{ek}=\sum _{j=1}^{N}p_{1j}{\frac {d_{k}}{d_{k}+d_{j}}}$ (2-114)

where N is the number of grain size classes.

Horizontal Sediment Mixing Coefficient

The horizontal sediment mixing coefficient, v_s, 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

$v_{s}=v_{t}/\sigma _{s}$(2-115)

where $\sigma _{s}$ is the Schmidt number and $v_{t}$ 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, $Q_{sed,tot}$ in kg/s, or as a fractional sediment transport rate,$Q_{sed,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 $Q{sed,k}=p_{1k}Q_{sed.tot}$ where $p_{1k}$ is the bed material fraction of the first layer at the boundary. The fractional sediment transport $q_{tk,B}$ at boundary cell B in kg/m/s is then calculated along the cell string according to

$q_{tk,B}={\frac {q_{tk,B}^{*}}{\sum \Delta l_{f}Q_{tk,N}^{*}}}Q_{sed,k}$(2-116)

where $q_{tk,B}^{*}$ is the potential sediment transport rate at boundary cell B and $\Delta l_{f}$ 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

${\frac {\partial (hC_{sal})}{\partial t}}+{\frac {\partial (hV_{j}C_{sal})}{\partial x_{j}}}={\frac {\partial }{\partial x_{f}}}\left(v_{sal}h{\frac {\partial C_{sal}}{\partial x_{j}}}\right)+S_{sal}$(2-117)

where

$C_{sal}$ = depth-averaged salinity [ppt or parts per thousand]

h = water depth [m]

$V_{j}$ = total flux velocity [m/s]

$v_{sal}$ = horizontal mixing coefficient $v_{sal}=v_{t}/\sigma _{s}\quad [m^{2}/s]$

$S_{sal}$ = 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 $U_{f}$ is zero between wet and dry cells and this eliminates advection transport. In addition, the salinity horizontal mixing coefficient at the cell face, ${\bar {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)

${\frac {\partial (2E_{sr}cw_{j})}{\partial x_{j}}}=-D_{sr}+f_{e}D_{br}$(2-118)

where

$E_{sr}$ = surface roller energy density [N/m]

$c$ = roller propagation speed [m/s]

$w_{j}=(cos\theta _{m},sin\theta _{m})$ is the wave unit vector [-]

$\theta _{m}$ = mean wave direction [deg]

$D_{sr}$ = roller energy dissipation [N/m/s]

$D_{br}$ = wave breaking dissipation (from the wave model) [N/m/s]

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

$D_{sr}={\frac {g2E_{sr}\beta _{D}}{c}}$(2-119)

where $\beta _{D}$ is the surface roller dissipation coefficient approximately equal to 0.05-0.1. The surface roller contribution to the wave stresses, $R_{ij}$ , is given by

$R_{ij}=2E_{sr}w_{i}w_{j}$(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).