TR-06-9

1Introduction to CMS-M2D

A coastal inlet connects an ocean, sea, or lake through a typically narrow landmass to the water body behind it, such as a bay, estuary, lagoon, or river. Hydrodynamic forcing contributing to the water exchange that maintains the coastal inlet channel may be the tide, river flow, wind, or seiching. A tidal inlet is the most common type of coastal inlet and is the waterway connecting the ocean and a bay, a lagoon, or a river entrance through which tidal and other currents, such as from rivers, flow. Longer period flows at coastal inlets, nearshore currents produced by waves, and wave action move sediment at many temporal and spatial scales, constrained by the sea bottom morphology, jetties, training structures, and navigation channel dredging. Change in morphology at coastal inlets in turn modifies the incident waves, circulation, and sediment pathways. Kraus (2006) gives an overview of morphologic response to coastal inlets and to engineering actions taken at inlets.

Advances in numerical modeling, increased speed of desktop computers, and reliable and powerful interfaces now make calculation of morphology change at coastal inlets feasible, despite the large number of interacting processes. Simulation of morphology change at inlets is an emerging area of applied research of interest to the U.S. Army Corps of Engineers (USACE) and its navigation mission at coastal inlets in design, operation, and maintenance of navigation channels and jetties, and in operating coastal inlets within the sediment-sharing system connecting the inlet, adjacent beaches, and bay.

Because of the wide range of problems to be addressed, both basic and applied, success in simulation of coastal inlet processes is promoted by exchange of information among interdisciplinary teams of coastal engineers and scientists working in government, private industry, and academia. A system approach is required that is based on a comprehensive numerical model with a code accessible to a wide range of researchers. At the same time, the model should be efficient in calculation so that it may be exercised in addressing practical problems. Experience with a variety of coastal inlets and engineering problems will call for continued improvement and advances in model capabilities.

The model M2D (Militello et al. 2004) is being developed under the Coastal Inlets Research Program (CIRP) conducted at the U.S. Army Engineer Research and Development Center (ERDC), Coastal and Hydraulics Laboratory (CHL) as a numerical simulation predictive tool with which to support multidisciplinary research teams and conduct practical projects at coastal inlets. CMS-M2D is a component of the CIRP’s Coastal Modeling System (CMS), formerly called the Inlet Modeling System, which is designed to calculate combined circulation (current and water-surface elevation), waves, and morphology change at inlets and nearby areas and operated on desktop computers through operation in the Surfacewater Modeling System (SMS) interface (Zundel 2000). In this interactive environment, CMS-M2D can simulate hydrodynamics, sediment transport, and morphology change at project level at coastal inlets and the connecting nearshore and bay. CMS-M2D is computationally efficient, easy to set up, and incorporates features required for many coastal engineering applications. CMS-M2D can be coupled to regional circulation models, such as the ADvanced CIRCulation (ADCIRC) model (Luettich et al. 1992), through boundary conditions (Militello and Zundel 2002a, 2002b), providing flexibility for larger scale applications and connectivity between models. At the same time, the code is accessible to CIRP researchers with expertise in areas other than computational fluid dynamics, as expected in interdisciplinary teams.

Prior to 2005, the model was called M2D. Because the model is central to the CMS in terms of its flexibility and functionality, it was renamed CMS-M2D in 2005. In this report, the model will be referred to as M2D or CMS-M2D, depending on if the reference is to the model prior to 2005 or to the present version.

The first report in the M2D series (Militello et al. 2004) documented Version 2.0, which calculated hydrodynamics under general forcing by tide, waves, and wind. The present report describes CMS-M2D Version 3.0, building upon the original technical documentation and describing sediment transport and morphology change calculations, other model advancements, changes to input and output formats, interface modifications, and example applications. An effort was made to provide a complete description of CMS-M2D in this report so that it is self-contained with respect to both the hydrodynamics and morphology change.

Scope of Model

CMS-M2D, Version 3.0, is a finite-volume numerical representation of the two-dimensional (2D) depth-integrated continuity and momentum equations of water motion. It also contains integrated representation of sediment (presently, sand-sized sediment) transport and morphology change through transport rate formulations, the advection-diffusion equation, and the sediment continuity equation for updating change in the sea bottom. Wave forcing is included in CMS-M2D through coupling with a wave model. In the CMS, radiation stresses from surface waves force currents and change the water-surface elevation and, if required, the calculated current can be input to the wave model to transform waves propagating on it. The governing equations, finite-volume approximations, representation of bottom and surface stresses, grid scheme, boundary conditions, other model features, and graphical interface are described in this report. Examples are given to demonstrate capabilities of CMS-M2D.

Calculation cells in CMS-M2D are defined on a staggered, rectilinear grid and can have constant or variable side lengths. Momentum equations are solved in a time-stepping manner first, followed by solution of the continuity equation, in which the updated velocities calculated by the momentum equations are applied. Optional sediment transport calculations are conducted on updated velocity and water-surface elevation values, together with wave properties if they are included in a specific simulation.

For support of engineering applications, features of CMS-M2D include flooding and drying, wind-speed dependent (time-varying) wind-drag coefficient, variably spaced bottom-friction coefficient, wave-stress forcing that can vary in space and time, efficient grid storage in memory, hot-start options, and the convenience of independently turning on or off the advective terms, mixing terms, wall friction, and sediment transport and morphology change calculations. Features of the sediment transport component include three transport formulation options, control over timing of transport calculations and morphology change calculations, representation of hard or non-erodible bottom, and representation of avalanching of sea-bottom slopes. Additionally, model output is written in formats convenient for importing into commercially available plotting packages. CMS-M2D operates exclusively in SI units and requires that input be provided in metric units.

CMS-M2D has been designed as a project-scale model that can be easily and quickly applied to examine engineering issues at coastal inlets such as navigation channel infilling, natural sand bypassing, consequences of mining ebb- or flood-tidal shoals, and changes in hydrodynamics and sediment transport in response to modifications to jetties. The model has been developed to maximize flexibility in grid specifications and forcing. For example, if forcing by water-surface elevation time series is prescribed as a boundary condition, the time interval between values in the time series does not have to be constant. Thus, if values were being supplied from a water-level gauge and the sampling rate was modified, the measurements could be prescribed as forcing without sub-sampling the data set.

Implementation within Surface-Water Modeling System

A graphical interface for CMS-M2D is implemented within the SMS, Versions 8.1 and higher (Zundel 2000). Features of the CMS-M2D interface cover grid development, control file specification, boundary condition and wind specification, coupling with regional and wave models, model runs, post-processing of results, and input and output visualization. The SMS provides flexibility in grid generation and modification by incorporating model-specific tools to assign bathymetric data sets for interpolating to the CMS-M2D grid, manually modifying depths and friction coefficients, adjusting cell sizes, and inserting or deleting calculation cells. Variably spaced grids can be readily developed with SMS tools that spatially scale the cell sizes over user-defined extents. The SMS contains dialogs for specification of boundary conditions. For some boundary types, forcing information can be input manually from within the SMS, or predefined information can be read from files. A model control dialog provides the user with a convenient interface for specifying timing control, model options, wind and wave forcing, sediment transport and morphology change options, and output options.

CMS-M2D can be driven by larger-domain circulation models, such as ADCIRC, through boundary specification capabilities contained within the SMS. The boundary conditions dialog allows access to solutions from larger-domain models that can be extracted and mapped to CMS-M2D boundaries. The user can decide whether to specify water-surface elevation or a combination of water-surface elevation and velocity as boundary input for CMS-M2D. This capability provides CMS-M2D with boundary conditions that preserve tidal phase and other variations calculated by a larger or regional scale model.

The SMS also allows automated coupling of CMS-M2D with the STeady state spectral WAVE model STWAVE (Smith et al. 1999, 2001) and the Wave Action Balance Equation with Diffraction (WABED) model (Mase 2001; Mase and Kitano 2000; Mase et al. 2005) through a Steering Module dialog, which is convenient for projects that require wave-stress forcing for CMS-M2D. The Steering Module allows the user to control the model coupling, giving flexibility in conducting simulations of wave-driven currents and the wave-current interaction.

Sediment transport and morphology change options are specified within the SMS model control dialog for CMS-M2D. This dialog allows the user to select from three sediment transport formulations and enter formulation-specific parameters. Timing controls for sediment transport rate calculations and for morphology change are included in the dialog.

Studies Performed with CMS-M2D

CMS-M2D has been developed and tested in basic and applied studies covering a wide range of bathymetric configurations, forcing, and project goals. Summaries of selected applications are described here to familiarize users with model capabilities and to establish model history.

M2D was first applied to investigate hydrodynamic responses to proposed engineering projects on the Texas coast. Inlet and bay hydrodynamics on this coast are frequently dominated by wind forcing, making this region a more complex study area than most inlet and bay systems on the Atlantic and Pacific Ocean coasts of the United States. Changes in circulation patterns in the Upper Laguna Madre and exchange of water with the Gulf of Mexico were calculated for proposed causeway modifications designed to improve evacuation safety (Brown, Militello, and Kraus 1995a, 1995b, 1995c) for the Kennedy Causeway that connects North Padre Island with the mainland. As a result of the analysis, the Kennedy Causeway was elevated to provide improved circulation to the Laguna Madre and a more reliable evacuation route from North Padre Island.

In an investigation of a proposed new inlet called the Southwest Corner Cut that would connect East Matagorda Bay to the Colorado River Navigation Channel, TX, M2D was applied to calculate water velocity changes in the navigation channel and Gulf Intracoastal Waterway and to determine whether the new inlet would promote closure of Mitchell’s Cut (Kraus and Militello 1996, 1999; Militello and Kraus 1998) an inlet and flood relief channel located on the opposite end of East Matagorda Bay.

In a basic study of the hydrodynamics of Baffin Bay, TX, M2D was applied to calculate the wind-driven circulation and set up and set down in this nontidal bay (Militello 1998, 2000). Model calculations also determined the contributions of nonlinear terms in the governing equations to harmonics of the oscillatory wind (sea breeze) forcing in Baffin Bay (Militello 1998; Militello and Kraus 2001).

In a channel infilling study for Bay Center, a harbor located inside Willapa Bay, WA (Kraus, Arden, and Simpson 2000), M2D calculated current velocity over a 1-month time interval from which bed elevation changes were computed (Militello 2002; Militello et al. 2003). The Bay Center application also demonstrated successful coupling of M2D with ADCIRC. This success, combined with CIRP goals, led to the development of automated M2D and ADCIRC (or other larger-domain model) coupling capability within the SMS (Militello and Zundel 2002a, 2002b).

As part of the CIRP’s efforts to assure model reliability, M2D calculated tide and wave-driven currents at an idealized inlet and ebb shoal for fair weather and storm waves (Militello and Kraus 2003). Interactions between wave-driven currents and the ebb jet varied, depending on the relative strength of the two current sources. For storm waves, calculated alongshore velocities were stronger and the surf zone was wider and farther from shore, as compared to the fair-weather waves. This study demonstrated coupling of M2D and a surface wave model (STWAVE), described in this volume and by Militello and Zundel (2003), and the capability to calculate sediment transport and morphology change that has since been advanced.

Recent applications include calculation of hydrodynamics and morphology change under combined tide and wave forcing at Shinnecock Inlet, NY (Bounaiuto and Militello 2004) and at the mouth of the Colorado River, TX (Lin, Kraus, and Barcak 2004). The study at Shinnecock Inlet compared measured change in bathymetry at the inlet and ebb shoal with calculations representing the same 6-month time frame (13 August 1997 to 28 May 1998). Observed major changes in morphology were shown to be represented by the model, including westward migration of the channel thalweg, scour of the inlet throat, and accretion of the ebb shoal. For the Colorado River study, six structural alternatives were evaluated for performance relating to sediment management at the inlet (river mouth) entrance. For the existing condition, M2D was demonstrated to reproduce sediment shoaling at the entrance and filling of a deposition basin located adjacent to a weir jetty.

Report Contents

This report describes Version 3.0 of CMS-M2D and gives instruction on development of a CMS-M2D project within the SMS. Chapter 2 covers the governing equations for water motion that are solved in the model. Chapter 3 similarly covers the governing equations for sediment motion and morphology change, as well as methods describing morphologic constraints that have been implemented in CMS-M2D. Chapter 4 describes the numerical approximations for water motion and defines the cell numbering method, and Chapter 5 describes the numerical approximations for sediment motion and morphology change. Chapter 6 reviews the types of boundary conditions available in CMS-M2D. Chapter 7 describes key model features, and Chapter 8 reviews the SMS interface for CMS-M2D and provides guidance for project development. Chapter 9 provides information on model setup and file structure. Chapter 10 provides example applications of CMS-M2D.

Appendices supplement the main text and contain more detailed and auxiliary information. Appendix A gives example input and output files. Appendix B provides references for CMS-M2D development and applications. Appendices C through F describe results of numerical model integrity tests. Appendix G demonstrates the wave-adjusted boundary condition, which modifies prescribed boundary condition values for the presence of wave forcing. Appendix H compares calculated and measured longshore currents from laboratory experiments. Appendix I demonstrates the capability for representation of hard bottom (non-erodible bottom).

2Governing Equations for Water Motion

CMS-M2D solves the 2D, depth-integrated continuity and momentum equations by applying a finite-volume method. Velocity components are calculated in two horizontal dimensions. This chapter introduces the governing equations and describes the respective terms and variables. Wind and bottom stress formulations are presented, as well as implementation of wave stresses.

The 2D depth-integrated continuity and momentum equations solved by CMS-M2D are:

[[Image:]](1)

[[Image:]](2)

[[Image:]](3)

where

h=still-water depth relative to a specific vertical datum

=deviation of the water-surface elevation from the still-water level

t=time

q_x=flow per unit width parallel to the x-axis

q_y =flow per unit width parallel to the y-axis

u=depth-averaged current velocity parallel to the x-axis

v=depth-averaged current velocity parallel to the y-axis

g=acceleration due to gravity

D_x=diffusion coefficient for the x direction

D_y=diffusion coefficient for the y direction

f=Coriolis parameter

<math>{{\tau }_{bx}}</math>

=bottom stress parallel to the x-axis

<math>{{\tau }_{by}}</math>

=bottom stress parallel to the y-axis

<math>{{\tau }_{wx}}</math>

=surface stress parallel to the x-axis

<math>{{\tau }_{wy}}</math>

=surface stress parallel to the y-axis

<math>{{\tau }_{Sx}}</math>

=wave stress parallel to the x-axis

<math>{{\tau }_{Sy}}</math>

=wave stress parallel to the y-axis

Component velocities are related to the flow rate per unit width as:

[[Image:]](4)

[[Image:]](5)

For the situation without waves, bottom stresses are given by:

[[Image:]](6)

[[Image:]](7)

where U = total current speed and C_b = empirical bottom-stress coefficient. The total current speed is:

[[Image:]](8)

The bottom-friction coefficient is calculated by:

[[Image:]](9)

where C = Chezy coefficient given by:

[[Image:]](10)

in which R = hydraulic radius, and n = Manning roughness coefficient.

In the presence of waves, the bottom stress contains contributions from both the quasi-steady current (as from tide, wind, and surface waves) and the bottom orbital motion of the waves. Assessment of an average bottom stress over the period of a surface wave, which requires numerical integration, would need to be performed at every grid point at every time step, a computation-intensive operation. Instead, here a square wave approximation for the wave is applied, which allows analytic estimation of the time average. It is given by the following approximation for the time-averaged bottom stress under combined currents and waves (Nishimura 1988):

[[Image:]](11)

[[Image:]](12)

in which [[Image:]] = wave angle relative to the x-axis, and [[Image:]] and [[Image:]] are given by:

[[Image:]](13)

[[Image:]](14)

where [[Image:]] = wave angular frequency, H = wave height, and k = wave number.

Surface wind stresses are given by:

[[Image:]](15)

[[Image:]](16)

where

[[Image:]]=wind drag coefficient

<math>{{\rho }_{a}}</math>

=density of air

<math>{{\rho }_{w}}</math>

=density of water

W=wind speed

 =wind direction

The convention for wind direction is specified to be 0 deg for wind from the east with angle increasing counterclockwise.

Wave stresses are calculated from spatial gradients in radiation stresses as:

[[Image:]](17)

[[Image:]](18)

where S_xx, S_xy, and S_yy are wave-driven radiation stresses. Radiation-stress tensor calculations are based on linear wave theory and computed within STWAVE, WABED, or other wave model for supplying information to CMS-M2D. They represent the summation of standard tensor formulations across the defined spectrum. For a coordinate system with the x-axis oriented normal to the shoreline, the tensor components are (Smith et al. 2001):

[[Image:]](19)

[[Image:]](20)

[[Image:]](21)

where

S_xx=flux of shore-normal momentum

S_xy=of shore-parallel momentum, or the shear component of the radiation stress

S_yy=flux of momentum alongshore

E=wave energy density

Wave models such as STWAVE typically calculate in coordinate systems such that the x-axis is normal to the shoreline (positive x directed toward shore) and the y-axis is parallel to the shoreline. For CMS-M2D applications, the radiation-stress tensor is rotated from the STWAVE grid and mapped onto the CMS-M2D grid by the SMS.

The Coriolis parameter is given by:

[[Image:]](22)

where [[Image:]] = angular frequency of the Earth's rotation, and [[Image:]] = latitude.

The depth-mean horizontal eddy viscosity coefficient D depends on the strength of mixing in the water column, which is a function of the acting processes at a particular location. If waves do not contribute significantly to mixing, D can be calculated as a function of the total water depth, current speed, and bottom roughness as (Falconer 1980):

[[Image:]](23)

where the subscript o denotes oceanic mixing. The form of the eddy viscosity coefficient given by Equation 23 produces a mixing term that is weakly nonlinear.

In the surf zone, waves contribute significantly to lateral mixing, and the eddy viscosity coefficient is expected to be a function of the wave properties. In CMS-M2D, surf zone mixing is represented by:

[[Image:]](24)

where [[Image:]] describes the lateral mixing below trough level (Smith, Larson, and Kraus 1993) and is expressed as (Kraus and Larson 1991):

[[Image:]](25)

where [[Image:]] = empirical coefficient representing the lateral mixing strength, and u_m = amplitude of the horizontal component of the wave orbital velocity at the bottom given by:

[[Image:]](26)

where T = wave period.

Mixing is a vigorous three-dimensional (3D) process in the surf zone not readily represented in a depth-averaged model. Such vertical processes cannot be directly described in a depth-average model, and the mixing coefficient in Equation 25 should be considered as a surrogate. Streaming from undertow advects turbulent eddies or mixing offshore, and this mixing can support a strong current seaward of the wave breakpoint (Putrevu and Svendsen 1992). With further distance seaward, the vertical current decreases as the volume of water increases, weakening the mixing and the longshore current.

The transition between surf zone mixing and oceanic mixing seaward of the breakpoint is represented in CMS-M2D by a weighted mixing coefficient specified as:

<math>D=\left( 1-{{\theta }_{m}} \right){{D}_{o}}+{{\theta }_{m}}{{D}_{w}}</math>(27)

where the weighting parameter <math>{{\theta }_{m}}</math> is given by:

<math>{{\theta }_{m}}={{\left( \frac{H}{h+\eta } \right)}^{3}}</math>(28)

The method of weighting and the weighting parameter are ad hoc and must be validated by comparison to data. Intuitively, a cubic dependence on water depth introduces the fluid volume.

To illustrate the control on the cross-shore distribution of the longshore current by wave friction and representation of mixing, four simulations were conducted that correspond to Visser (1982) Case 4. The Visser laboratory data and CMS-M2D grid and parameter values are described in Appendix H with one difference being that the Manning’s n value for velocities discussed in this section was 0.013. The four simulations consisted of combinations of wave friction specification and mixing coefficient. Simulations denoted as “No Wave Friction” did not include waves as a source of friction, whereas “Wave Friction” denoted the bottom stress as defined by Equations 11 and 12. Mixing coefficients were calculated as “Wave Mixing,” which indicated application of Equation 24, and “Weighted Mixing,” which denoted the combined oceanic and surf zone mixing described by Equation 27.

Calculated and measured longshore current for Visser Case 4 are compared in Figure 1. Values of distance increase toward the offshore direction. Peak calculated current speed varies between 0.37 and 0.40 m/s and measured peak speed is 0.39 m/s. Thus, combinations of both wave friction and representation of the mixing coefficient exert control on the maximum current speed. Greatest variation among the calculated speeds is located between the peak current and the weakening current with distance offshore. Calculation with the weighted mixing coefficient and wave friction produced a curve in which the magnitude and slope compare well with measurements. All remaining calculations show deviation in magnitude and curvature from measured values in the offshore tail.

Accurate estimation of the wind-stress coefficient is necessary for calculating wind-induced flow in shallow water (depths [[Image:]]5 m). A short review describing the development of the variable wind-stress coefficient implemented in CMS-M2D is presented here.

[[Image:]]

Figure 1. Comparison of longshore current speed values for Visser Case 4

The flux of momentum into the water column from the wind is specified as a function of the wind speed and wind-stress coefficient. Anemometer height above the water or land surface must be considered for accurate estimation of the wind-stress coefficient because wind speed varies logarithmically in the vertical direction for unstratified atmospheric conditions. A general law for the vertical wind profile is given by (Charnock 1955; Hsu 1988):

[[Image:]](29)

where

W_z=wind speed at height Z above the surface

Z₀=surface roughness

[[Image:]]=friction velocity

[[Image:]]=von Kármán constant

The friction velocity is defined in terms of the surface wind stress [[Image:]] as (Hsu 1988):

[[Image:]](30)

where K_m = eddy viscosity coefficient, and z is the vertical dimension.

Under the assumption of nearly neutral atmospheric stability, the wind-stress coefficient over water at a height 10 m above the water surface, C₁₀, can be expressed as (Hsu 1988):

[[Image:]](31)

where W₁₀ = wind speed in meters per second at 10 m above the surface. Equation 31 includes the assumption of fully developed seas. For measurements of wind speed made at heights other than 10 m, an approximation for the 10-m wind speed is (Shore Protection Manual 1984):

[[Image:]](32)

The surface wind-stress Equations 15 and 16, as implemented in CMS-M2D, apply the 10-m wind-drag coefficient and wind speed to compute the surface stress. Thus, values of C₁₀ and W₁₀ defined by Equations 31 and 32 are specified as C_d and W in Equations 15 and 16.

Equations 31 and 32 have been successfully applied in modeling studies investigating motion in strong wind-forced shallow embayments with typical depths ranging from 1 to 4 m (Brown, Militello, and Kraus 1995a; Kraus and Militello 1999; Militello 2000). Calculated water level and current compared well with measurements, indicating that these formulations are satisfactory for estimation of the vertical wind profile and wind stress in shallow-water systems.

3Governing Equations for Sediment Motion and Methods for Morphologic Constraints

As a user-specified option, CMS-M2D calculates sediment transport rates and resultant changes in water depth (bottom elevation) through gradients in the transport rates. In CMS-M2D Version 3, three transport rate formulations are available:

1. Watanabe (1987) total load formulation.
2. Lund-CIRP (Camenen and Larson 2005, 2006) total load formulation (combined suspended and bed load).
3. Advection-diffusion (AD) transport for suspended load coupled with either the van Rijn (1998) or the Lund-CIRP formulation for bed load. The AD equation applies the reference concentration and sediment diffusivity from either the van Rijn or Lund-CIRP formulation as a user option.

The sediment transport formulas are applicable with or without the presence of waves. Change in bottom elevation is calculated by the sediment continuity equation in which transport rates computed from a specified formulation are entered. The Lund-CIRP formulas have the properties of smooth transitions between areas of breaking and non-breaking waves and between areas of transport by current only and by a current and waves.

CMS-M2D accounts for the two morphologic constraints of: (a) hard bottom (non-erodible substrate) and (b) bottom avalanching if the slope exceeds a specified value. The approach taken to treat hard bottom allows for sand to accumulate over the hard bottom if depositional conditions are present, and it also allows material overlying the hard substrate to erode while preserving sediment volume. Avalanching is invoked if the bottom slope between two cells exceeds a specified critical slope. The avalanching algorithm applies an iterative approach to move material down slope until the critical slope is no longer exceeded, while conserving sediment volume.

The transport rate formulas that can be applied in CMS-M2D Version 3 are based on the shear stress concept, implying that bed roughness and resulting friction factors are key parameters to estimate when computing transport rates. A sediment transport formula is often developed (i.e., calibrated against data) based on specific equations for calculating bed roughness and associated shear stress. Thus, when calculating the transport rate in CMS-M2D with the three implemented formulas, slightly different approaches are taken to determine the shear stress. In the following sections, the methods are discussed for obtaining the shear stress for each of the sediment transport formulas. Also, it is noted that the roughness and associated shear stresses within the sediment transport calculations are different from those in the CMS-M2D hydrodynamic calculations. The computed current field and properties of surface waves calculated from a separate model are input to the roughness and shear stress calculations for the sediment transport, but there is presently no feedback of roughness or shear stress from the sediment transport component of the model to the hydrodynamic or wave calculations.

CMS-M2D calculates the current velocity in two horizontal coordinate directions (u, v). If there is negligible contribution from the waves to the transporting velocity (e.g., sinusoidal waves), the net transport is in the direction of the resulting current vector and has magnitude of

<math>{{U}_{c}}=\sqrt{{{u}^{2}}+{{v}^{2}}}</math>

. Thus, waves will contribute to the mobilization and stirring of sediment, but a mean current is needed to produce a net transport. The Watanabe total load formulation and the AD equation compute sediment transport rates along the x and y axes. The Lund-CIRP total load formula computes transport in the direction of the current. In order to obtain the Lund-CIRP calculated transport in the two horizontal coordinate directions, the transport vector is split up into x- and y-directed components.

The following sections describe the sediment transport equations applied in CMS-M2D, followed by the methodologies implemented to treat hard bottom and avalanching.

Watanabe Formulation

The total load sediment transport rate of Watanabe (1987) is given by:

<math>{{q}_{tot}}=A\left[ \frac{\left( {{\tau }_{b,\max }}-{{\tau }_{cr}} \right){{U}_{c}}}{{{\rho }_{w}}g} \right]</math>(33)

where

<math>{{q}_{tot}}</math>=total load (both suspended and bed load)

<math>{{\tau }_{b,\max }}</math>

=maximum shear stress at the bed

[[Image:]]=shear stress at incipient sediment motion

<math>{{U}_{c}}</math>=depth averaged current velocity

[[Image:]]=density of water

[[Image:]]=acceleration of gravity

<math>A</math>=empirical coefficient typically ranging from 0.1 to 2

The shear stress at incipient motion is:

[[Image:]](34)

where

[[Image:]]=sediment density

[[Image:]]=median grain diameter

[[Image:]]=critical Shields parameter

For the case of currents only (no waves) the bed shear stress

<math>{{\tau }_{c}}</math>

is:

<math>{{\tau }_{c}}=\frac{1}{8}{{\rho }_{w}}{{f}_{c}}{{U}_{c}}^{2}</math>

(35)

where [[Image:]] is the current friction factor and

<math>{{\tau }_{b,\max }}={{\tau }_{c}}</math>

in Equation 33. The current friction factor is given as (van Rijn 1998):

<math>{{f}_{c}}=0.24{{\log }^{-2}}\left( 12\frac{d}{{{k}_{sd}}} \right)</math>(36)

in which d = total depth defined as

<math>d=h+\eta </math>

, and k_sd = Nikuradse equivalent sand grain roughness obtained from:

<math>{{k}_{sd}}=2.5{{d}_{50}}</math>(37)

If waves are present, the maximum bed shear stress [[Image:]] is given by a data-based method of (Soulsby 1997):

[[Image:]](38)

where [[Image:]] = mean shear stress by waves and current over a wave cycle,[[Image:]] = wave bed shear stress, and [[Image:]] = angle between the waves and the current. The mean wave and current shear stress is:

<math>{{\tau }_{m}}={{\tau }_{c}}\left[ 1+1.2{{\left( \frac{{{\tau }_{w}}}{{{\tau }_{c}}+{{\tau }_{w}}} \right)}^{3.2}} \right]</math>

(39)

The wave bed shear stress is given by:

<math>{{\tau }_{w}}=\frac{1}{2}{{\rho }_{w}}{{f}_{w}}U_{w}^{2}</math>

(40)

where [[Image:]] = wave friction factor, and

<math>{{U}_{w}}</math>

= wave orbital velocity amplitude. The wave friction factor is (Nielsen 1992):

[[Image:]](41)

where R is the relative roughness defined as:

<math>R=\frac{{{A}_{w}}}{{{k}_{sd}}}</math>(42)

where A_w is the semi-orbital excursion defined as:

<math>{{A}_{w}}=\frac{{{U}_{w}}T}{2\pi }</math>(43)

Transport rates are set to zero at the beginning of the simulation, then updated at each sediment transport time step.

Lund-CIRP Formulation

The Lund-CIRP formulation is implemented into CMS-M2D in two modes. One mode is calculation of total load through combining bed load and suspended load transport. The second mode is implementation in the AD transport equation. Bed roughness and friction factors applied in the Lund-CIRP formulation are first introduced, followed by descriptions of the bed load and suspended load formulations, respectively.

Bed roughness and friction factors

Bottom roughness k_s is assumed to consist of three components, namely grain related roughness k_sd, form-drag roughness k_sf, and sediment-related roughness k_ss (Soulsby 1997). Total roughness is obtained as the linear sum of these components:

<math>{{k}_{s}}={{k}_{sd}}+{{k}_{sf}}+{{k}_{ss}}</math>

(44)

where the grain-related roughness is given by Equation 37.

Form-drag roughness is primarily a function of the ripple properties and calculated as (Soulsby 1997):

<math>{{k}_{sf}}=7.5\frac{H_{r}^{2}}{{{L}_{r}}}</math>

(45)

in which H_r = ripple height, and L_r = ripple length, which are computed from:

<math>\begin{align}

& {{L}_{r}}=1000{{d}_{50}} \\

& {{H}_{r}}=\frac{{{L}_{r}}}{7} \\

\end{align}</math>

(46)

for a current (Soulsby 1997) and from:

[[Image:]](47)

for waves (van Rijn 1993), where

<math>{{\psi }_{w}}</math>

= mobility parameter that is defined by:

<math>{{\psi }_{w}}=\frac{U_{w}^{2}}{(s-1)g{{d}_{50}}}</math>

(48)

and s = specific density given by:

<math>s=\frac{{{\rho }_{s}}}{{{\rho }_{w}}}</math>

(49)

The sediment-related roughness is given by Wilson (1966, 1989):

<math>{{k}_{ss}}=5{{d}_{50}}{{\theta }_{i}}</math>

(50)

where

<math>{{\theta }_{i}}</math>

= Shields parameter for current or waves (i = c, w, respectively, for current and waves). Expressions for the Shields parameter for currents

<math>{{\theta }_{c}}</math>

and for waves

<math>{{\theta }_{w}}</math>

are given by:

<math>{{\theta }_{c}}=\frac{{{\tau }_{c}}}{\rho \left( s-1 \right)g{{d}_{50}}}</math>

(51)

<math>{{\theta }_{w}}=\frac{{{\tau }_{w}}}{\rho \left( s-1 \right)g{{d}_{50}}}</math>

(52)

Equation 50 has to be solved simultaneously with the expression for the shear stress, because the roughness depends on the stress. Explicit predictive equations were developed using approximating polynomials to the exact solution avoid time-consuming iteration in calculating the shear stress.

Based on the roughness estimate, friction factors for current f_c and waves f_w are separately computed according to Soulsby (1997) and Swart (1974), respectively as:

<math>{{f}_{c}}=2{{\left( \frac{\kappa }{1+\ln \left( {{k}_{s}}/30d \right)} \right)}^{2}}</math>

(53)

[[Image:]](54)

where

<math>\kappa </math>

= Von Karman constant (taken to be 0.4). Shear stress from the current is given by:

<math>{{\tau }_{c}}=\frac{1}{2}{{\rho }_{w}}{{f}_{c}}U_{c}^{2}</math>

(55)

whereas the maximum shear stress from the waves is estimated from:

<math>{{\tau }_{w}}=\frac{1}{2}{{\rho }_{w}}{{f}_{w}}U_{w}^{2}</math>

(56)

and the mean shear stress is obtained from:

<math>{{\tau }_{wm}}=0.5\,{{\tau }_{w}}</math>

(57)

where a sinusoidal wave is assumed. The corresponding shear velocities for current and waves, respectively, are defined as:

<math>{{u}_{*c}}=\sqrt{\frac{{{\tau }_{c}}}{{{\rho }_{w}}}}</math>

(58)

for the current, and:

<math>{{u}_{*w}}=\sqrt{\frac{{{\tau }_{w}}}{{{\rho }_{w}}}}</math>

(59)

for the waves (based on the maximum shear stress).

Bed load

A general transport formula for bed load

<math>{{q}_{b}}</math>

under combined waves and current was developed by Camenen and Larson (2005):

<math>\begin{align}

& \frac{{{q}_{bw}}}{\sqrt{(s-1)gd_{50}^{3}}}={{a}_{w}}\sqrt{{{\theta }_{net}}}{{\theta }_{cw,m}}\exp \left( -b\frac{{{\theta }_{cr}}}{{{\theta }_{cw}}} \right) \\

& \frac{{{q}_{bn}}}{\sqrt{(s-1)gd_{50}^{3}}}={{a}_{n}}\sqrt{{{\theta }_{cn}}}{{\theta }_{cw,m}}\exp \left( -b\frac{{{\theta }_{cr}}}{{{\theta }_{cw}}} \right) \\

\end{align}</math>

(60)

where subscripts w and n refer to the wave direction and the direction normal to the waves, respectively, a and b are coefficients, and

<math>{{\theta }_{cw,m}}</math>

and

<math>{{\theta }_{cw}}</math>

are the mean and maximum Shields parameter for waves and current combined, respectively, not including the bed form roughness. The coefficient a_n = 12 for the transport component normal to the wave direction (current-driven transport). The coefficient b describes initiation of motion and is specified as b = 4.5.

The quantities

<math>{{\theta }_{net}}</math>

and

<math>{{\theta }_{cn}}</math>

represent the net contribution of the shear stress to the transporting velocity during a wave cycle in the direction parallel and normal to the waves, respectively, defined as:

<math>{{\theta }_{net}}={{\theta }_{cw,on}}+{{\theta }_{cw,off}}</math>

(61)

in which:

<math>\begin{align}

& {{\theta }_{cw,on}}=\frac{1}{{{T}_{wc}}}\int\limits_{0}^{{{T}_{wc}}}{\frac{1}{2}\frac{{{f}_{cw}}{{\left( {{u}_{w}}(t)+{{U}_{c}}\cos \phi \right)}^{2}}}{(s-1)g{{d}_{50}}}dt} \\

& {{\theta }_{cw,off}}=\frac{1}{{{T}_{wt}}}\int\limits_{{{T}_{wc}}}^{T}{\frac{1}{2}\frac{{{f}_{cw}}{{\left( {{u}_{w}}(t)+{{U}_{c}}\cos \phi \right)}^{2}}}{(s-1)g{{d}_{50}}}dt} \\

\end{align}</math>

(62)

and

<math>{{\theta }_{cn}}=\frac{1}{2}{{f}_{c}}\frac{{{\left( {{U}_{c}}\sin \phi \right)}^{2}}}{(s-1)g{{d}_{50}}}</math>

(63)

where T_wc = period of positive flow, T_wt = period of negative flow, and u_w = time-varying horizontal bottom orbital velocity (

<math>T={{T}_{wc}}+{{T}_{wt}}</math>

). Analytical solutions to these equations are available for simple expressions of u_w. The friction factor for combined wave and current flow is calculated from Madsen and Grant (1976):

<math>{{f}_{cw}}={{X}_{v}}{{f}_{c}}+\left( 1-{{X}_{v}} \right){{f}_{w}}</math>

(64)

where

<math>{{X}_{v}}=\left| {{U}_{c}} \right|/\left( \left| {{U}_{c}} \right|+{{U}_{w}} \right)</math>

and f_c and f_w are the friction coefficients for current only and waves only, respectively, not including the bed form roughness. The Shields parameter that determines the amount of sediment mobilized for transport is obtained as:

<math>{{\theta }_{cw,m}}={{\left( \theta _{c}^{2}+\theta _{w,m}^{2}+2{{\theta }_{w,m}}{{\theta }_{c}}\cos \phi \right)}^{1/2}}</math>

(65)

and the Shields parameter applied in the initiation of motion term is:

<math>{{\theta }_{cw}}={{\left( \theta _{c}^{2}+\theta _{w}^{2}+2{{\theta }_{w}}{{\theta }_{c}}\cos \phi \right)}^{1/2}}</math>

(66)

where

<math>{{\theta }_{c}}</math>

and

<math>{{\theta }_{w}}</math>

are the Shields parameters for current only and waves only, respectively, and the index m denotes the mean value. For sinusoidal waves,

<math>{{\theta }_{w,m}}={{\theta }_{w}}/2</math>

.

Values of the transport rate coefficients were obtained through calibration against a large data set from laboratory and field experiments involving a mean current, oscillatory waves (sinusoidal and asymmetric), and combined waves and current. A wide range of current and wave velocities, grain sizes, and roughness conditions were covered in this combined data set. The following relationship was developed for the transport coefficient in the wave direction (Camenen and Larson 2005):

<math>\begin{align}

& {{a}_{w}}=6+6{{X}_{t}} \\

& {{X}_{t}}=\frac{{{\theta }_{c}}}{{{\theta }_{c}}+{{\theta }_{w}}} \\

\end{align}</math>

(67)

In many applications, the waves are assumed to be sinusoidal, having no asymmetry. Thus, the contribution to the transporting velocity from the waves is negligible, implying that only the current moves the material. Waves are still important for mobilization of the sediment. For such conditions, Equation 60 simplifies to:

<math>\frac{{{q}_{bc}}}{\sqrt{(s-1)gd_{50}^{3}}}={{a}_{n}}\sqrt{{{\theta }_{c}}}{{\theta }_{cw,m}}\exp \left( -b\frac{{{\theta }_{cr}}}{{{\theta }_{cw}}} \right)</math>

(68)

where the transport q_bc is obtained in the direction of the current, and the transport normal to the current is zero. (There is no velocity component to move the material in this direction.)

Suspended load

The suspended load q_s is calculated based on the assumptions of an exponential concentration profile and a constant velocity over the water column to yield (Camenen and Larson 2006):

<math>{{q}_{s}}={{U}_{c}}{{c}_{R}}\frac{\varepsilon }{{{w}_{s}}}\left( 1-\exp \left( -\frac{{{w}_{f}}d}{\varepsilon } \right) \right)</math>

(69)

where w_f = sediment fall speed, c_R = reference concentration, and

<math>\varepsilon </math>

= sediment diffusivity (or mixing coefficient). The transport q_s is taken to be in the direction of the current because the waves are assumed to produce a zero net drift and not contribute to the transport. In the calculation of c_R and

<math>\varepsilon </math>

, the total shear stress should be applied, including the bed form roughness. The reference concentration is given by:

<math>{{c}_{R}}={{A}_{cR}}{{\theta }_{cw,m}}\exp \left( -b\frac{{{\theta }_{cr}}}{{{\theta }_{cw}}} \right)</math>

(70)

where the coefficient <math>{{A}_{cR}}</math> is determined by the following relationship:

<math>{{A}_{cR}}=3.5\cdot {{10}^{-3}}\exp \left( -0.3{{d}_{*}} \right)</math>

(71)

and the dimensionless grain size is defined by

<math>{{d}_{*}}={{\left( (s-1)g/{{\nu }^{2}} \right)}^{1/3}}{{d}_{50}}</math>

. The sediment diffusivity is related to the energy dissipation according a generalization of the treatment of Battjes (1975):

<math>\varepsilon ={{\left( \frac{{{D}_{e}}}{\rho } \right)}^{1/3}}d</math>

(72)

where D_e is a total effective dissipation given by:

<math>{{D}_{e}}=k_{c}^{3}{{D}_{ec}}+k_{w}^{3}{{D}_{ew}}+k_{b}^{3}{{D}_{eb}}</math>

(73)

in which k_c, k_w, and k_b are coefficients. The dissipations from bottom friction due to current D_ec and from bottom friction due to waves D_ew are expressed as:

<math>{{D}_{ec}}={{\tau }_{c}}{{u}_{c*}}</math>

(74)

<math>{{D}_{ew}}={{\tau }_{w}}{{u}_{w*}}</math>

(75)

and the dissipation due to wave breaking D_eb is directly provided from a wave transformation model (e.g., STWAVE. WABED). The coefficient k_b corresponds to an efficiency factor and was found to be 0.017, whereas k_c and k_w are related to the Schmidt number that expresses the ratio between diffusivity of sediment and water (i.e., momentum). These two coefficients are given by:

<math>{{k}_{i}}=\frac{1}{6}{{\sigma }_{i}}\kappa </math>

(76)

where

<math>{{\sigma }_{i}}</math>

= Schmidt number and the subscript i = c or w denotes current and waves, respectively. The Schmidt number is calculated from the empirical relationships:

[[Image:]](77)

where A₁ = 0.7 and A₂ = 3.6 for current, and A₁ = 0.09 and A₂ = 1.4 for waves. In the case of coexistent waves and current, a weighted Schmidt number is applied according to:

<math>{{\sigma }_{cw}}=X_{v}^{5}{{\sigma }_{c}}+\left( 1-X_{v}^{5} \right){{\sigma }_{w}}</math>

(78)

for both waves and current.

Bed load and suspended load transport rates are set to zero at the beginning of the simulation, then updated at each sediment transport time step.

Advection-Diffusion (AD) Equation

Calculations of suspended load and bed load are conducted separately through solution of the AD equation for suspended load. This method calculates the suspended load in the water column through the lateral and vertical movement of the sediment.

Basic equations

Sediment transport consists of two components, bed load and suspended load. The bed load and suspended load transport rate formulas are based on the local shear stress. However, in the situation where the suspended sediment concentration rapidly changes in time and space, such as at a river mouth, inlet entrance, navigation channel, and in the vicinity of structures, the suspended load cannot be determined solely by the local forcing conditions. In this case, suspended sediment concentration is governed by the AD equation.

The AD equation implemented in CMS-M2D is obtained from the continuity of depth-averaged suspended sediment transport (Figure 2) as:

<math>\frac{\partial (Cd)}{\partial t}+\frac{\partial (C{{q}_{x}})}{\partial x}+\frac{\partial (C{{q}_{y}})}{\partial y}=\frac{\partial }{\partial x}\left( {{K}_{x}}d\frac{\partial C}{\partial x} \right)+\frac{\partial }{\partial y}\left( {{K}_{y}}d\frac{\partial C}{\partial y} \right)+P-D</math>

(79)

in which

C=depth-averaged sediment concentration

d=total depth (= h+)

h=still-water depth

=deviation of the water-surface elevation from the still-water level

t=time

qx=flow per unit width parallel to the x-axis (= ud)

qy=flow per unit width parallel to the y-axis (= vd)

u=depth-averaged current velocity parallel to the x-axis

v=depth-averaged current velocity parallel to the y-axis

Kx=sediment diffusion coefficient for the x direction

Ky=sediment diffusion coefficient for the y direction

P=sediment pick-up rate (upward sediment flux)

D=sediment deposition rate (downward sediment flux)

The suspended sediment concentration can be solved from Equation 79 at a sediment transport time step.

[[Image:]]

Figure 2. Definition of continuity (conservation of volume) of sediment transport

Bed level change is obtained from the sediment continuity equation:

<math>\frac{\partial h}{\partial t}=\frac{1}{1-p}\left( \frac{\partial {{q}_{bx}}}{\partial x}+\frac{\partial {{q}_{by}}}{\partial y}+P-D \right)</math>

(80)

where q_bx = bed load transport rate parallel to the x-axis, q_by = bed load transport rate parallel to the y-axis, and p = porosity of sediment.

In Equation 80, the terms

<math>{\partial {{q}_{bx}}}/{\partial x}\;+{\partial {{q}_{by}}}/{\partial y}\;</math>

include the bed load component, and the term

<math>P-D</math>

corresponds to the suspended load component. Figure 2 shows the definition of sediment continuity in a water column, where q_sx and q_sy are the suspended transport rate in the x and y directions, respectively, and a is the height of the bed load layer. The transport rates q_sx and q_sy may be defined as follows from Equation 79:

<math>{{q}_{sx}}=C{{q}_{x}}-{{K}_{x}}d\frac{\partial C}{\partial x}</math>

(81)

<math>{{q}_{sy}}=C{{q}_{y}}-{{K}_{y}}d\frac{\partial C}{\partial y}</math>

(82)

Sediment concentration and bed load and suspended load transport rates are set to zero at the beginning of the simulation, then updated at each sediment transport time step.

Pick-up and deposition rates

The AD model calculates the bed level change due to suspended load from the difference between pick-up rate and deposition rate in Equation 80. The pick-up rate and the deposition rate are also applied as the bottom boundary condition in Equation 79. The boundary conditions are specified at an arbitrary level a above the mean bed level:

<math>P={{\left. -\varepsilon \frac{\partial c}{\partial z} \right|}_{z=a}}={{c}_{a}}{{w}_{f}}</math>

(83)

<math>D={{c}_{0}}{{w}_{f}}</math>

(84)

where c = equilibrium concentration of suspended sediment at a given elevation, and z = vertical coordinate. Both c_a and c₀ are reference concentrations defined at z = a. Because the upward flux of sediment depends on the bed shear stress, c_a is determined from the bed shear stress calculated from the local hydrodynamic conditions. Representation of c_a within CMS-M2D is dependent on selection of either the van Rijn or Lund-CIRP models. The downward sediment flux depends on the concentration in the upper water column; therefore, c₀ is specified from solution of Equation 79.

Assuming that the suspended concentration is in equilibrium (

<math>{\partial c}/{\partial t}\;=0</math>

,

<math>{\partial u}/{\partial x}\;={\partial v}/{\partial y}\;=0</math>

,

<math>{\partial c}/{\partial x}\;={\partial c}/{\partial y}\;=0</math>

), then the basic equation for suspended sediment concentration can be written:

<math>c{{w}_{f}}+\varepsilon \frac{\partial c}{\partial z}=0</math>

(85)

This equation can be solved analytically by applying an appropriate mixing coefficient , and the vertical profile of suspended concentration is obtained in the following form:

<math>c(z)={{c}_{0}}F(z)</math>

(86)

where F(z) is a function of the vertical concentration distribution. The relationship between the reference concentration c₀ and the depth-averaged concentration C is:

<math>C=\frac{1}{d-a}\int_{a}^{d}{{{c}_{0}}F(z)}dz</math>

(87)

Thus, c₀ may be written in the following form by introducing a conversion function _d:

<math>{{c}_{0}}=\frac{C}{\frac{1}{d-a}\int_{a}^{d}{F(z)}dz}=\frac{C}{{{\beta }_{d}}}</math>

(88)

The bed level change due to suspended load is based upon the difference of the two types of reference concentration:

<math>P-D=({{c}_{a}}-{{c}_{0}}){{w}_{f}}</math>

(89)

This equation implies that erosion occurs if c_a > c₀, and accretion occurs if c_a < c₀.

In the AD model, three methods of specifying c_a and c₀ (that is, _d) are implemented. Two of the methods are based on the van Rijn formula (van Rijn 1985) and one on the Lund-CIRP formula (Camenen and Larson 2006). Table 1 summarizes general features of the methods.

van Rijn formula

The reference concentration c_a in the van Rijn formula is given by:

<math>{{c}_{a}}=0.015\frac{{{d}_{50}}}{a}{{\left( \frac{{{\tau }_{s,\max }}-{{\tau }_{cr}}}{{{\tau }_{cr}}} \right)}^{1.5}}d_{*}^{-0.3}</math>

(90)

where a = reference height and _s,max = maximum bed shear stress given by Equation 38.

Table 1Features of Calculation of Pick-up and Deposition Rates
Method	Reference Concentration	Conversion Parameter	Comments
Exponential Profile	Eq. 90	Eq. 104 and 105	Fast computation. Tends to overestimate sediment transport rate. Can be used for some tests.
Van Rijn Profile	Eq. 90	Solves Eq. 85 numerically with Eq. 103 (Runge-Kutta 4th)	Requires substantial computing time. Provides the same results as van Rijn (1985).
Lund-CIRP Profile	Eq. 70	Eq. 104 and 72	Fast computation. Newly developed sediment transport formula.

The current-related shear stress is calculated from:

<math>{{\tau }_{c}}=\frac{{{\rho }_{w}}{{\kappa }^{2}}}{{{\left[ 1+\ln \left( {{{{{k}'}}_{s}}}/{30d}\; \right) \right]}^{2}}}U_{c}^{2}</math>

(91)

The wave-related shear stress is given as:

<math>{{\tau }_{w}}={{\rho }_{w}}\frac{{{f}_{w}}}{2}U_{w}^{2}</math>

(92)

<math>{{f}_{w}}=\exp \left[ 5.5{{\left( {{{A}_{w}}}/{{{{{k}'}}_{s}}}\; \right)}^{-0.2}}-6.3 \right]</math>

(93)

where <math>k_{s}^{'}</math> = roughness height defined as:

<math>{{{k}'}_{s}}={{k}_{sd}}+{{k}_{ss}}</math>

(94)

in which k_sd and k_ss are calculated from Equations 37 and 50, respectively.

Bed concentration in the van Rijn model is defined at the height a as:

<math>a=\max \left( 0.5{{H}_{r}}\,,\,0.01d \right)</math>

(95)

where H_r = ripple height. If ripples are present, the total roughness height is modified as:

<math>{{k}_{s}}={{{k}'}_{s}}+{{k}_{sf}}</math>

(96)

where, k_sf = 20(H_r²/L_r) is specified for the van Rijn formula, and k_sf = 7.5(H_r²/L_r) is specified for the Lund-CIRP formula. The ripple dimension is obtained by selecting the largest ripple height for the case of current or waves (Equations 46 and 47).

The conversion parameter _d to determine c₀ is obtained from the vertical mixing coefficient. Van Rijn (1985) proposed a distribution of the mixing coefficients for only current or waves according to Figure 3. The current-related mixing coefficient is given by:

[[Image:]](97)

where u_*c = current-related bed shear velocity expressed as:

<math>{{u}_{*c}}=\frac{\kappa }{\left[ -1+\ln \left( \frac{30d}{{{k}_{s}}} \right) \right]}{{U}_{c}}</math>

(98)

and <math>\vartheta </math> is a coefficient obtained from:

[[Image:]](99)

[[Image:]]

Figure 3. Vertical distributions of mixing coefficient due to current and waves

The wave-related mixing coefficient is:

[[Image:]](100)

where

<math>\begin{align}

& {{\varepsilon }_{w,\text{bed}}}=0.00065{{D}_{*}}{{\alpha }_{br}}\delta {{u}_{w}}\stackrel{\scriptscriptstyle\to}{\leftarrow} \\

& \\

& {{\varepsilon }_{w,\max }}=0.035{{\alpha }_{br}}\frac{hH}{T} \\

\end{align}</math>

(101)

[[Image:]](102)

in which the parameter

<math>\delta </math>

Z = height from the bed given by

<math>\delta =3{{H}_{r}}</math>

, H = significant wave height, and T = significant wave period. If waves and current coexist, the combined mixing coefficient is given by:

<math>{{\varepsilon }_{cw}}=\sqrt{\varepsilon _{c}^{2}+\varepsilon _{w}^{2}}</math>

(103)

By applying the expression for _cw, the concentration profile can be derived from Equation 85, and the conversion parameter _d calculated. In the AD-model based on the van Rijn equations, two methods are implemented to calculate _d. One is based on an exponential profile employing a depth-averaged mixing coefficient

<math>{{\bar{\varepsilon }}_{cw}}</math>

, and the other uses the original van Rijn profile, where the _d is obtained by numerical integration of Equation 85. Assuming an exponential profile of suspended sediment concentration, _d is obtained analytically as:

<math>{{\beta }_{d}}=\frac{1}{d-a}\frac{{{{\bar{\varepsilon }}}_{cw}}}{{{w}_{f}}}\left[ 1-\exp \left( -\frac{{{w}_{f}}}{{{{\bar{\varepsilon }}}_{cw}}}(d-a) \right) \right]</math>

(104)

where

<math>{{\bar{\varepsilon }}_{cw}}=\frac{1}{d-a}\int_{a}^{d}{{{\varepsilon }_{cw}}dz}</math>

(105)

Lund-CIRP formula

Reference concentration and sediment diffusivity calculated by the Lund-CIRP formula may also be applied in the AD-model. The formulas used for c_R and ε are presented in the section describing the Lund-CIRP formula.

Horizontal diffusion coefficient

Horizontal diffusion coefficients K_x and K_y for the depth-averaged sediment concentration in the AD equation are calculated by the equation proposed by Elder (1959):

<math>\begin{align}

& {{K}_{xo}}=5.93{{u}_{*c}}d \\

& {{K}_{yo}}=5.93{{u}_{*c}}d \\

\end{align}</math>

(106)

In the surf zone, waves contribute significantly to lateral mixing, and the eddy viscosity coefficient is expected to be a function of the wave properties. For the horizontal diffusion coefficients for AD equation, surf zone mixing proposed by Kraus and Larson (1991) is applied as well as for hydrodynamic calculation.

<math>{{K}_{w}}={{\varepsilon }_{L}}</math>

(107)

where [[Image:]] describes the lateral mixing below trough level (Smith, Larson, and Kraus, 1993) and is expressed as (Kraus and Larson 1991):

<math>{{\varepsilon }_{L}}=\Lambda {{u}_{w}}H</math>

(108)

where [[Image:]] = empirical coefficient representing the lateral mixing strength.

The transition between surf zone diffusion and oceanic diffusion seaward of the breakpoint is represented in CMS-M2D by a weighted mixing coefficient specified as:

<math>\begin{align}

& {{K}_{x}}=(1-{{\theta }_{m}}){{K}_{xo}}+{{\theta }_{m}}{{K}_{w}} \\

& {{K}_{y}}=(1-{{\theta }_{m}}){{K}_{yo}}+{{\theta }_{m}}{{K}_{w}} \\

\end{align}</math>

(109)

where the weighting parameter <math>{{\theta }_{m}}</math> is given by:

<math>{{\theta }_{m}}={{\left( \frac{H}{h+\eta } \right)}^{3}}</math>

(110)

Bed load transport rate

Bed load transport rate may be calculated with the van Rijn (1998) or Lund-CIRP (Camenen and Larson 2005) formula. In the van Rijn formula, bed load transport rate is obtained from, not including bed roughness:

<math>{{q}_{b}}=0.1{{d}_{50}}{{u}_{*c}}{{\left( \frac{{{\tau }_{s,\max }}-{{\tau }_{cr}}}{{{\tau }_{cr}}} \right)}^{1.5}}d_{*}^{-0.3}</math>

(111)

The Lund-CIRP formula for bed load is described previously in this chapter.

Scaling Factors for Bed Load and Suspended Load Transport Rates

Adjustments to the bed load and suspended load transport rates can be assigned for both the Lund-CIRP total load formulation and the AD equation. Scaling factors are applied directly to transport rate values as coefficients independently for the bed load transport rate and for the suspended load transport rate. Thus, the bed load transport rate can be scaled separately from the suspended load transport rate. Recommended ranges for the scaling factors are 0.1 to 4. A scaling factor of 1 will not modify the transport rate, and this value is the default. Values less than 1 will reduce the transport rate, and values greater than 1 will increase it. Scaling factors are sometimes found necessary because of uncertainties in forcing and initial conditions, allowing adjustment of model predictions of sediment transport rates and morphology change to observations.

Sediment Continuity Equation

Depth change for the total load formulations and for bed load transport from the AD equation is calculated by the sediment continuity equation, in which time-averaged transport rates are applied. This approach requires that two time intervals be defined, one over which to compute instantaneous transport rates dt_sed and the other over which to compute averaged transport rates and morphology change dt_morph. Directional instantaneous transport rates

<math>{{q}_{tot,x}}</math>

and

<math>{{q}_{tot,y}}</math>

are calculated at intervals of dt_sed and averaged over the time interval dt_morph to obtain the time-averaged transport rates

<math>{{\bar{q}}_{tot,x}}</math>

and

<math>{{\bar{q}}_{tot,y}}</math>

for the x and y directions, respectively. Enhancement of down-slope transport and reduction of up-slope transport rates is accomplished by calculation of modified time-averaged transport rates

<math>\bar{q}_{tot,x}^{'}</math>

and

<math>\bar{q}_{tot,y}^{'}</math>

as (Watanabe 1987):

<math>\bar{q}_{tot,x}^{'}={{\bar{q}}_{tot,x}}+{{D}_{s}}\left| {{{\bar{q}}}_{tot}} \right|\frac{\partial h}{\partial x}</math>(112)

<math>\bar{q}_{tot,y}^{'}={{\bar{q}}_{tot,y}}+{{D}_{s}}\left| {{{\bar{q}}}_{tot}} \right|\frac{\partial h}{\partial y}</math>(113)

in which [[Image:]] = empirical slope coefficient with typical range of 5 to 30, and

<math>{{\bar{q}}_{tot}}</math>

is the time-averaged transport rate. These formulations add an effective diffusion of transport to the continuity equation for bed change.

The sediment continuity equation is solved to compute depth change at the time interval dt_morph and for total load formulations is given by:

<math>\frac{dh}{dt}=\left( \frac{1}{1-p} \right)\left( \frac{\partial \bar{q}_{tot,x}^{'}}{\partial x}+\frac{\partial \bar{q}_{tot,y}^{'}}{\partial y} \right)</math>(114)

At forcing boundaries, the gradient in transport rate across the boundary cell is set to zero.

For the AD-model, the sediment continuity equation is provided by Equation 80. The time intervals to compute instantaneous transport rates and to compute averaged transport rates and morphology change are consistent with those applied in the total load model, that is, dt_sed and dt_morph respectively. Also, the bed load components in Equation 80 are modified by the slope effect as:

<math>{{\bar{{q}'}}_{bx}}={{\bar{q}}_{bx}}+{{D}_{s}}\left| {{{\bar{q}}}_{b}} \right|\frac{\partial h}{\partial x}</math>

(115)

<math>{{\bar{{q}'}}_{by}}={{\bar{q}}_{by}}+{{D}_{s}}\left| {{{\bar{q}}}_{b}} \right|\frac{\partial h}{\partial y}</math>

(116)

where

<math>{{\bar{q}}_{b}}=\sqrt{\bar{q}_{bx}^{2}+\bar{q}_{by}^{2}}</math>

. Finally, the sediment continuity equation is solved to compute depth change at the time interval dt_morph and is rewritten as:

<math>\frac{dh}{dt}=\frac{1}{1-p}\left( \frac{\partial {{{{\bar{q}}'}}_{bx}}}{\partial x}+\frac{\partial {{{{\bar{q}}'}}_{by}}}{\partial y}+\bar{P}-\bar{D} \right)</math>

(117)

which contains both bed load and suspended load components. At forcing boundaries, the gradient in transport rate across the boundary cell is set to zero.

Hard-Bottom Methodology

Representation of hard (non-erodible) substrate provides the capability to simulate mixed bottom types within a single simulation. The default bottom type in CMS-M2D is erodible, meaning that sediment is available to be transported from that area. Areas having hard bottom are specifically designated as such in the model, and the depth of the hard substrate is specified, allowing representation of either exposed or sediment-covered hard bottom. The one-dimensional (1D) hard-bottom calculation approach applied in the GENESIS shoreline change model (Hanson 1989; Hanson and Kraus 1989) to represent the control of seawalls on shoreline evolution (Hanson and Kraus 1985, 1986) was extended to two dimensions and implemented in CMS-M2D (Hanson and Militello 2005).

The presence of hard bottom exerts a constraint on water depth because the bottom cannot erode. Representation of hard bottom applies an adjustment to the potential transport rates (calculated by a transport formulation with no restriction) to account for the presence of hard substrate. Treatment of hard bottom was developed to satisfy the following four properties:

a.Areas having hard bottom cannot erode below the level of the hard substrate.

b.Sediment volume is conserved.

c.The direction of sediment transport at a location with hard bottom is the same as that of the direction of potential local transport.

d.Correction to transport rates is always done in the direction of potential sediment transport. If hard bottom is exposed, the reduction in actual transport as compared to the situation without hard bottom is always transmitted to neighboring areas in the direction of transport as governed by the hydrodynamic forcing.

The general procedure for treating hard bottom is:

a.Calculate sediment transport and morphology change without hard-bottom restriction.

b.Assess whether hard-bottom cells have eroded below the hard bottom depth. If so, apply a correction procedure that adjusts transport rates based on available sediment at hard-bottom cells and corrects depths so that the hard-bottom constraint is not violated. The correction procedure makes adjustments at all cells in which transport rates are controlled by hard bottom and can include adjustments at cells that do not have hard bottom.

The correction procedure involves four steps:

a.Identify all cells that have transport at all four faces directed out of the cell, called “Minus4” cells. These Minus4 cells, if they exist, must constitute the first stage of corrections, as they are only 'deliverers' of sediment to all adjacent cells, not 'recipients' of material.

b.If the hard bottom is exposed, all transport rates out of the Minus4 cells are reduced. A list of all adjacent cells, with or without hard bottom, that are influenced by exposed hard-bottom cells is generated for later correction to depths and transport rates.

c.Depths at all adjacent cells are then recalculated to account for the reduced influx of sediment. Cells with hard bottom must then be checked for a possible violation of the hard-bottom constraint and corrected as necessary. For each corrected cell, all sediment transport rates out of the cell are in turn corrected, and each of the next generation of affected cells adjacent to it is added to the list of cells to be corrected. Transport rates into a cell are never changed (as it is an outgoing transport from the adjacent cell that is corrected in that cell).

d.After all cells have been checked, the process is repeated to make sure that the hard-bottom constraint is not violated in any cell.

Corrections to transport rates and depths are based upon the potential transport rate and amount of sand available for erosion. If, after solution of the sediment continuity equation (Equation 114), the updated depth [[Image:]] is greater than the hard-bottom level [[Image:]], the outgoing transport rates are corrected to match the criterion [[Image:]]. Over one morphologic time step, the total amount of material<math>{{h}_{avail}}</math>, expressed as a vertical quantity that may leave the cell without violating the hard-bottom constraint is:

<math>{{h}_{avail}}=({{h}_{hb}}-{h}')+\frac{{{q}_{in}}\,dt}{dx\,dy}</math>

(118)

where q_in = total influx of sediment from the cell during the morphologic time step. The correction factor K is given as:

[[Image:]](119)

where q_out = total outflux of sediment from the cell during the morphologic time step. To ensure that <math>{h}'={{h}_{hb}}</math> ingoing transport rates are left unchanged, and outgoing transport rates are multiplied by the correction factor K.

A schematic example of the correction procedure is provided for a symmetric computational grid (Figure 4). The domain contains a 9-by-9 grid of calculation cells. Rows (parallel to the x-axis) are numbered 1 through 9, and columns (parallel to the y-axis) are indicated by A through I. Initial bathymetry level is specified to be <math>z=0</math> at all cells, and the hard-bottom level is at <math>z=5</math> (units) at all cells (z is positive downwards from the original bottom level at <math>z=0</math>).

Arrows across the cell walls indicate the initial (potential) transport rates. These are of three different magnitudes: 2 units/m/s (black arrows), 5 units/m/s (blue arrows), and 12 units/m/s (green arrows). For demonstration purposes we assume that the cell sizes in the x-direction and the y-direction are the same and that a transport rate of 1 unit/m/s during one time step corresponds to a bottom elevation change of 1.

Numbers inside each cell indicate the potential depth change (positive numbers mean erosion) in each cell, which is a measure of the magnitude of the depth change if the hard bottom were not there. Thus, in each cell with a depth change that is greater than 5, the hard bottom is exposed. As a result, the actual depth change is limited to 5. Once this depth is reached, the transport out of the cell is reduced owing to limited sediment availability.

In the model, this reduction is accomplished through a correction procedure. The following section explains how this correction is performed through three successive correction loops in a hierarchical and systematic order.

[[Image:]]

Figure 4. Schematic hard-bottom grid and potential transport rates

Correction loop 1:

a.Find cells with transport rates directed out of the cell through all four faces. These cells represent the ultimate starting points. Cell E5 is the only cell in this example having outward-directed transport rates through all four faces.

b.Calculate the new bottom level equaling 8 units, which means that the hard bottom at 5 units is exposed.

c.Reduce all transport rates out of cell E5 in proportion to the sediment deficit. The potential loss of sediment was 8 units while the available volume in cell E5 was 5 units. Thus, the outgoing transport rates will be reduced from 2 to 2 x 5/8 = 1.25 units/m/s.

d.Recalculate depths of all adjacent cells, D5, E4, F5, and E6, based on modified transport rates.

e.Hard bottom is now exposed in cells E4 and E6.

f.Reduce transport rates out of these two cells (E4 and E6).

g.Recalculate depths of cells adjacent to E4 and E6 that receive sediment. These cells are D4, E3, F4 and D6, F6, E7. No additional hard-bottom cells are exposed. Thus, the deficit associated with cell E5 does not spread further.

Correction loop 2:

h.Find start cells from where sand-deficit might spread by searching through the grid from one end to the other. Start cells are those that are located on a hard bottom, but do not receive sediment from any other hard-bottom cell.

i.Find exposed cells (defined as cells where the hard bottom is exposed) and identify by flag IEXP = 1. Start the search in cell A1. This cell is exposed, but sediment is coming in from B1 that is also exposed. The direction of influence will be in the direction of transport, i.e., from B1 to A1. Thus, A1 cannot be a start cell.

j.Next, proceed to cell A2. This cell is exposed but, as above, it is receiving sediment from cell A1, which is also exposed. Thus, cell A2 is not a start cell.

k.The same situation prevails for all cells in the A column. Thus, no cells in the A column are start cells.

l.Next, proceed to cell B1. This cell is exposed. It is only getting sand from B2, but that cell is not exposed, so no correction can come from there. Thus, cell B1 is a start cell.

m.All remaining cells are checked the same way, and all start cells are added as necessary to the list of cells to be corrected (here referred to as the correction list).

n.Correct the first cell in the correction list, cell B1. Add adjacent cells receiving sediment from B1 to the correction list, i.e., cells A1 and C1.

o.Proceed through all start cells, correct outgoing transport rates, and add to the correction list all cells receiving sediment from the start cells.

Correction loop 3:

p.The next step is to examine all cells connected to start cells. Transport rates need to be corrected if hard bottom is exposed. Adjacent cells receiving sediment are added to the correction list at the next level of correction cells.

Through this procedure, cells will always be corrected in a hierarchic order in which all cells of one level are corrected before proceeding to the next level. Corrections are always propagated in the direction of transport without ambiguity regarding which cell affects another. Figure 5 shows the final bathymetry where all (positive) depth changes are less than 5. The figure also shows the transport rates that were unchanged within this correction procedure (indicated by their original colors) together with those transport rates that were changed (indicated by red arrows). It should be noted, however, that the red arrows represent transport rates with very different magnitudes.

Through the correction loops, transport rates modified by the hard-bottom algorithm are not restricted to cells specified as having non-erodible substrate. Adjustments to transport rates take place at hard-bottom and erodable cells and the spatial extent of the modifications is dependent on the spatial influence of the hard bottom under the transport conditions.

[[Image:]]

Figure 5. Corrected hard-bottom depths and adjusted transport rates

Avalanching Methodology

Controls on bottom slope are implemented through an avalanching algorithm. After each calculation of morphologic change, every cell in the CMS-M2D grid is evaluated to determine if a critical bottom slope has been exceeded. If so, avalanching takes place to reduce the slope. By this method, realistic bathymetric slopes are maintained over the CMS-M2D domain. The method extends an algorithm developed for the Storm-Induced Beach (SBEACH) change model (Larson and Kraus 1989) to two dimensions.

The procedure for determining if avalanching is needed is entered after each calculation of morphologic change. Four slopes are computed for each cell in the grid using the updated depth values. The slopes are computed between each cell and its four neighbors (top, right, bottom, and left neighbors). Starting at the deepest cell in the grid, the bottom slopes are compared to the critical slope, which is the maximum slope allowed before avalanching takes place. If more than one cell has the same depth, cells are checked in order of cell identification number. If the critical slope is exceeded, avalanching takes place with the ending slope being equal to the residual slope after shearing. In the avalanching process, each avalanching calculation includes only the two neighboring cells that exceed the critical slope. Depths of avalanched cells are updated, and the slopes for each of the four neighboring cells recomputed. If more than one slope has exceeded the critical slope, the steepest is avalanched first as it is assumed that it would be the first to reach a critical state. If there is more than one slope that exceeds the critical limit by the same amount, then the slopes are avalanched from north to west in a clockwise fashion.

After an avalanche has occurred, the slopes are recomputed and checked to ensure they have not remained in a critical state. If a critical slope remains, the profile is again avalanched and the process repeated. After all critical slopes have been eliminated, the slopes of all affected cells are recomputed and the algorithm proceeds to the next deepest cell in the grid to repeat the process. All active cells are checked, whether they are wet or dry. After all cells have been checked and, if necessary, avalanched, slopes are again checked against the critical slope. If critical slopes remain, the algorithm finds the deepest cell and the entire process is repeated. If no profiles exceed the critical state, the avalanching process is complete.

A bias can result for cases where there are multiple cells with the same depth or multiple profiles that exceed the critical state by the same degree. In practice, however, it is unlikely that there will repeatedly be multiple cells with equivalent depths or slopes as the values are computed with a high degree of precision.