3Numerical Methods
Contents
 1 Numerical Methods
 1.1 Overview
 1.2 Computational Grid
 1.3 Spatial Discretization
 1.4 Temporal Discretization
 1.5 Cellface interpolation
 1.6 Cellface Gradient
 1.7 Cellcentered Gradient
 1.8 Reconstruction, Monotonicity, and Slope Limiters
 1.9 Advection Schemes
 1.10 Source/Sink Term
 1.11 Assembly of Algebraic Equations
 1.12 Implicit Relaxation
 1.13 Iterative Solvers
 1.14 Convergence and TimeStepping
 1.15 Ramp Function
 2 Hydrodynamics
 3 Salinity Transport
 4 Sediment Transport and Morphology Change
 4.1 Transport Equations
 4.2 Bed Change Equations
 4.3 Bed Material Sorting
 4.4 Avalanching
 4.5 Hard bottom
 4.6 Coupling of Sediment Transport, Bed Change, and Sorting Equations
 4.7 Summary of Sediment Transport SemiCoupling Procedure:=
 4.8 Spatial Interpolation and Extrapolation
 4.9 Temporal Interpolation and Prediction
Numerical Methods
Overview
CMSFlow has both implicit and explicit solution schemes. The explicit solver is designed for dynamic problems with extensive wetting and drying that require small computational time steps, while the implicit solver is intended for simulating tidal and waveinduced circulation at tidal inlets, navigation channels, and adjacent beaches. A detailed description of the numerical formulation of the explicit solver of CMSFlow can be found in Buttolph et al. (2006) and is not repeated here. The sections below specifically refer to the implicit solver of CMSFlow.
The implicit solver uses the SIMPLEC (SemiImplicit Method for Pressure Linked Equations Consistent) algorithm (Van Doormal and Raithby 1984) on a nonstaggered grid to handle the coupling of water level and velocity. Primary variables u, vvelocity, and water level are stored on the same set of grid points, and fluxes at cell faces are determined using a Rhie and Chow (1983) type momentum interpolation method (Wu et al. 2011). The explicit solver uses a staggered grid with velocities at the cell faces and the water levels and water depths at the cell centers (Buttolph et al. 2006a). CMSFlow also calculates salinity, sediment transport, and morphology change. The governing equations for hydrodynamics and sediment and salinity transport have similar forms which can be written as a general transport (advectiondiffusion) equation. In order to avoid redundant derivations of discretized equations, the discretization of the general transport equation is described in this chapter, and the same discretization may be applied to all of the transport equations. Then, the specific solution procedures for hydrodynamics, sediment transport, and bed change are introduced.
Computational Grid
The explicit time scheme in CMSFlow supports uniform and nonuniformly spaced Cartesian grids while the implicit version of CMSFlow supports generic Cartesian grids which can be uniform, nonuniform, or telescoping. The telescoping locally refines the mesh by splitting a cell into subcells. The only requirement imposed by the numerical methods is that the cells must have a rectangular shape. Additional requirements are imposed by the user interface which limits the variety of types of Cartesian grids to help simplify the grid generation and avoid grid quality issues. The following requirements are applied:
 Cells can only be subdivided into four subcells.
 Only two neighboring cells are allowed in the same direction (i.e., north, south, east, and west).
 Cells may have a maximum of six neighbors.
 Refinement levels must be spaced by at least one cell apart (i.e., cells that share the same corner must be one refinement level apart).
Examples of violations of the above four requirements are shown in sequential order in Figure 31 from a to d. Limiting the number of neighboring cells to six avoids excessive cell refinement and limits the band width of the coefficient matrix for the linearlized system of equations solved by the implicit solution scheme. The last two requirements avoid having excessive cell refinement which can cause numerical instabilities. The first requirement simplifies the grid generation process but may be relaxed in future versions. The last requirement is enforced for grid quality purposes and a more smoothly varying grid resolution.
The CMSFlow grids generated in the SMS interface can be classified as uniform, nonuniform Cartesian grids, regular telescoping, and stretched telescoping grids (Figure 32).
One important aspect of incompressible flow models is the location of primary variables: velocity and pressure (water level). On a staggered grid, the pressure (water level) is located at the center of cells and the u and vvelocities are located along the faces of cells (Harlow and Welsh 1965; Patankar 1980). On a nonstaggered grid, all of the primary variables are located at the center of cells. A staggered grid can more easily eliminate the checkerboard oscillations when compared to a nonstaggered grid; however, a nonstaggered grid involves a simpler source code and can minimize the number of coefficients that must be computed and stored during a simulation because many of the terms in the equations are approximately equal. In particular, a staggered grid is more complicated when handling the interface between coarse and fine cells where five or sixface control volumes are used. Therefore, a nonstaggered (collocated) grid approach is adopted for CMSFlow, with a Rhie and Chow (1983) momentum interpolation technique used to eliminate the checkerboard oscillations. Figure 33 shows the location of primary variables and the 5 and 7point stencils (computational molecule).
The data structure for a grid can be approached in three ways: 1) blockstructured, 2) hierarchical tree, and 3) unstructured. The blockstructured approach divides the domain into multiple blocks, and each block is treated as structured. The block structured approach requires a special treatment between blocks to ensure mass and momentum balance using this approach. The hierarchical tree approach is memory intensive and requires parentchild relationships and a tree traverse to determine mesh connectivity. For the unstructured approach, all cells are numbered in a 1D sequence, and tables are used to determine the connectivity of neighboring cells. Among these three approaches, the unstructured approach is the most simple and is therefore applied in CMSFlow. Computational cells are numbered in an unstructured manner via a 1D index array. Inactive cells (permanently dry) are not included in the 1D index array to save memory and computational time. All active computational cells are numbered sequentially. It is noted that using an unstructured data approach does not limit the model from using matrix solvers designed for structured grids. For convenience with handling boundary conditions, each boundary cell has a neighboring ghost cell outside of the computational domain.
General Transport Equation=
In order to avoid redundant derivations of discretized equations, discretization of a general transport equation is outlined and described below. All of the governing equations are some form of this transport equation; therefore, the same discretization may be applied to all of the governing equations. The general transport equation is given by

(1) 
where:
 = general scalar
 h = total water depth [m]
 = correction factor []
 = total flux velocity [m/s]
 = diffusion coefficient for []
 S = source/sink term which includes all remaining terms.
Note that in the case of the continuity and momentum equations, is equal to 1 and respectively. In the case of sediment and salinity transport, is equal to and , respectively. The correction factor (β) is equal to the totalload correction factor and equal to 1 for all other equations.
Spatial Discretization
A controlvolume technique is used in which the governing equations are integrated over a control volume to obtain an algebraic equation that can be solved numerically. Integration of Equation (1) over a control volume (shaded areas in Figure 33) yields the following:

(2) 

(3) 

(4) 
where
 = outward unit vector normal to cell face f[]
 = outward cell face velocity [m/s]
 = value of at the cell centroid
 = total water depth at the cell face centroid [m]
 = interpolated total water depth at cell face f [m]
 = outward normal gradient of at cell face f
 = interpolated diffusion coefficient for
 = area of cell P [m/s]
 = length of cell face f [m]
 = source/sink term.
In the above equations, the Gauss Divergence Theorem has been used to convert the area integral to a boundary integral. In addition, the area ingegrals have been evaluated assuming variable vary linearly within cells. The cell face velocity, , is calculated using a momentum interpolation method similar to that of Rhie and Chow (1983) and is described in Hydrodynamics of the present Chapter.
Temporal Discretization
The general transport equation is rewritten as
(5) 
where includes all the remaining terms. For stability and efficiency, a fully implicit timestepping scheme is used in the form
(6) 
where is a weighting factor between 0 and 1. For , the scheme reduces to the firstorder backward Euler scheme, and with , the scheme reduces to the secondorder backward scheme (Ferziger and Peric 1997). The superscripts indicate the time step levels, with n+1 being the next timestep and n being the previous timestep.
Cellface interpolation
The general formula for estimating the cellface value is given by
(7) 
where is a linear interpolation factor, is the gradient operator in the direction parallel to face f, and is the distance from the cell center to the ghost point (O) parallel to the cell face (f) (see Figure 34).
The subscripts and indicate the components parallel and normal to the cell face. By definition, . Note that for neighboring cells without any refinement, is equal to zero and the above equation is consistent with nonrefined cell faces. The linear interpolation factor is defined as

(8) 
where is the coordinate of f perpendicular to the face and is the cell dimension perpendicular to the face f.
Cellface Gradient
A linearly exact secondorder approximation for the normal gradient at cell face f is calculated using the auxiliary node concept of Ferziger and Peric (1997):
(9) 
where the subscripts P and N refer to two neighboring cells, is the distance between cells P and N, normal to the cell face (Figure 34), and is the gradient operator in the direction parallel to face f. Ham et al. (2002) compared the auxiliary node formulation to the fullyunstructured discretization proposed by Zwart et al. (1998) for the viscous terms and found that the auxiliary node formulation is significantly more stable.
Cellcentered Gradient
The cellcentered gradient operator is calculated using the Gauss’ Divergence Theorem as
(10) 
The equation above is second order and conservative for regular and nonuniform grids.
Reconstruction, Monotonicity, and Slope Limiters
Variables are linearly reconstructed within the cells as
(11) 
where is the cellaverage value specified at the cell centroid, is the distance vector from the cell centroid to any location within the cell, and is the cellcentered gradient. The reconstruction is second order and conservative in the sense that . The linear reconstruction is used when interpolating cellface values (Equation 7) and calculating cellface gradients (Equation 13). If the reconstruction satisfies the local maximum principle
(12) 
then no new extrema are created within the cell, and the solution is monotonic. Figure 35 shows two examples of linear reconstruction with and without slope limiters to ensure monotonicity.
For nontelescoping grids, the slope limiter is calculated as
(13) 
where R is the ratio between two consecutive slopes
(14) 
Here, the secondorder van Leer (1979) limiter is used because of its smoothness. Figure 36 shows a comparison of three different common limiters. The slope limiter is applied in each direction separately.
For unstructured grids the slope limiters described above are difficult to implement because of the difficulty in defining forward and backward differences. For telescoping grids the following Limited Central Difference (LCD) slope limiting procedure of Hubbard (1999) is applied:
(15) 
where . In the procedure outlined by Hubbard (1999), a scalar limiter is then calculated as . Here, a directional limiter is applied as , which is less dissipative. Finally, the cellcentered gradient is limited as
(16) 
where is the unlimited gradient.
Advection Schemes
Hybrid Scheme
The hybrid scheme is composed of a firstorder upwind scheme and a secondorder central difference scheme. The cellface advective value is given by
(17) 
where the subscripts D and C indicate the downstream and upstream values, and is the Peclet number at the cell face in which ,. In the hybrid scheme, when the Peclet number is larger than 2, the firstorder upwind value is used; otherwise, the secondorder central difference value is used.
Exponential Scheme
The exponential scheme interpolates the face value using an exact solution to the 1D steady advectiondiffusion equation:
(18) 
The exponential scheme has automatic upwinding and is stable but is usually less than second order.
Hybrid Linear/Parabolic Scheme
The Hybrid Linear/Parabolic Approximation (HLPA) scheme of Zhu (1991) may be written as
(19) 
where the subscripts D, C, and U indicate the first downstream and first and second upstream cells, respectively. The normalized variable ) is determined based on the formulation of Jasak et al. (1999)
(20) 
where . The HLPA scheme is second order. Choi et al. (1995) found that the HLPA scheme has similar accuracy to the thirdorder SMARTER (Sharp and Monotonic Algorithm for Realistic Transport Efficiently Revised) and LPPA (Linear and PiecewiseParabolic Approximation) schemes but is simpler and more efficient (Shin and Choi 1992; Choi et al. 1995).
Source/Sink Term
The source/sink term is linearized as (Patankar 1980)
(21) 
The coefficient is required to be nonpositive for stability.
Assembly of Algebraic Equations
Assembly refers to the process of combining terms to create a linear system of algebraic equations. The algebraic equation for each cell is obtained by first combining or assembling all of the terms. Then, the continuity equation is multiplied by and is subtracted from the transport equation. The resulting discretized equation for cell P is
(22) 
where the subscript N refers to the neighboring cell sharing cell face, and are linear coefficients for and . The last term, (b), contains all the remaining terms. Applying a similar equation for all of the internal cells of a grid results in a system of algebraic equations. This set of equations are referred to as the discretized governing equations.
Implicit Relaxation
The righthand side of the algebraic system of equations (b) generally contains deferred corrections and source terms which must be computed iteratively. This iteration loop is referred to as the outer loop, and the loop within the linear iterative solver (see Section Iterative Solvers) is referred to as the inner loop. Underrelaxation is often applied to stabilize the convergence of the outer iteration loop. Underrelaxation may be applied as in which is an underrelaxation parameter . The superscript m indicates the previous iteration, and superscript m+1 indicates the new relaxed value. A better approach is used here to introduce the underrelaxation directly in the algebraic system of equations (Majumdar 1988; Ferziger and Peric 1997) as
(23) 
This is referred to implicit underrelaxation and has the effect of making the coefficient matrix more diagonally dominant, thus improving convergence.
Iterative Solvers
The selection of an iterative solver is one key issue concerning the overall performance of the model. CMS has three iteration solvers available and they are described in more detail below: 1) GMRES variation, 2) BiCGStab, and 3) GaussSeidel. The default iterative solver for CMS is a variation of the GMRES (Generalized Minimum RESidual) method (Saad 1993) and is used to solve the algebraic equations. The original GMRES method (Saad and Schultz 1986) utilizes the Arnoldi process to reduce the coefficient matrix to the Hessenburg form and minimizes the norm of the residual vector over a Krylov subspace at each iterative step. The variation of the GMRES method used here allows changes in preconditing at every iteration step Saad (1993). The Incomplete Lower Upper Factorization ILUT (Saad 1994) is used as the preconditioner to speedup convergence. The GMRES solver is applicable to symmetric and nonsymmetric matrices and leads to the smallest residual for a fixed number of iterations. However, the memory requirements and computational costs become increasingly expensive for larger systems.
The BiCGStab (BiConjugate Gradient Stabilized) iterative solver is also a Krylov subspace solver and is applicable to symmetric and nonsymmetrix matrices (Saad 1996). BiCGStab also uses ILUT as a preconditioner (Saad 1994). The BiCGStab method can be viewed as a combination of the standard Biconjugate Gradient solver where each iterative step is followed by a restarted GMRES iterative step. One advantage of the BiCGStab iterative solver is that the memory requirements are constant for each iteration and there are less computational costs when compared to the GMRES method (for less than four iterations).
The SIP (Strongly Implicit Procedure) iterative solver uses an Incomplete Lower Upper decomposition, which approximates the exact Lower Upper decomposition (Stone 1968). The method is specifically designed for algebraic systems of equations derived from partial differential equations. The implementation here is for a 5point stencil and, therefore, only applies to nontelescoping grids.
The ICCG (Incomplete Cholesky preconditioned Conjugate Gradient) iterative solver is applicable to symmetric matrices such as the pressure correction equation (Ferziger and Peric 2002). The implementation here is also for a 5point stencil and, therefore, can only be applied to nontelescoping grids
The simplest iterative solvers implemented here are the pointimplicit GaussSeidel solvers with or without SuccessiveOverRelaxation. The SuccesiveOverRelaxation may speed up convergence but can also lead to model divergence (Patankar 1980). Even though the GaussSeidel method requires more iterations for convergence, the overall efficiency may be higher than the GMRES and BiCGStab because each iteration is computationally inexpensive, and the code is efficiently parallelized. However, based on experience and testing, the GMRES and BiCGStab methods are usually more robust and perform better for large timesteps.
Convergence and TimeStepping
During the iterative solution process, error is calculated and used to determine if a solution has converged, diverged, or stalled at an error below a set tolerance threshold. An estimate of the error in solving the general algebraic equation is given by
(24) 
The norm of the residual errors is given by
(25) 
Since this value depends on the total number of cells, the final statistic (referred to as the residual) that is used for estimating the model convergence is obtained by dividing the norm by the square root of the number of cells:
(26) 
where is the number of computational cells. Simply stated is referred to as the normalized residual error and the superscript refers to the iteration number. Also is calculated for each variable that is solved at each iteration step of the solution process. Each equation has default maximum tolerances for determining if the solution has converged, diverged, or stalled. The maximum number of iterations that is imposed is set equal to M. A minimum of 5 iterations are required for the hydrodynamic equations, and a minimum of M/2 iterations are required for the sediment transport equations. Table 31 lists the default criteria to determine whether the iterative solution procedure has converged, requires a reduced time step or has diverged.
Table 31. Default criteria to determine whether the iterative solution procedure has converged, diverged, or requires a reduced time step.
Current velocity,m/s  If R^{m} < 1x10^{7}  If R^{M} > 1.0x10^{3}  If R^{M} > 1.0x10^{2} or U_{i} > 10.0 
Pressure correction, m^{2}/s^{2} 
If R^{M} < 1x10^{8}  If R^{M} > 1.0x10^{4}  If R^{M}>1.0x10^{3} or P > 50.0 
Totalload concentration, kg/m^{3} 
If R^{M}< 1x10^{8}  None  If R^{M} > 1.0x10^{3} or Ctk < 0 
Salinity, ppt  If R^{m}< 1x10^{6}  None  If Csal < 0 
For the implicit model, the time steps for the hydrodynamics, sediment and salinity transport are the same in order to avoid mass conservation problems and for simplicity. If any of the timestep reduction criteria are met, then the timestep is reduced by half and a minimum number of two timesteps are calculated at the newly reduced time step. If the last timestep converged properly, then the timestep is increased. The maximum time step allowed is equal to the userspecified initial timestep.
Ramp Function
For most coastal applications, the model is initialized from a cold start, which means that the water level and current velocities are initially set to zero. When the model is initialized with anything other than zeros, this is referred to as a hot start. The ramp period allows the model to slowly transition from the specified initial conditions without shocking the system. The ramp function is defined here as
(27) 
where t is the simulation time and is the ramp duration. The ramp function provides a smooth function for transitioning from the initial condition and is plotted in Figure 37.
Figure 3.7. Ramp function used in CMS
The ramp function is applied to the model forcing conditions, including the wave forcing, surface wind, sediment concentration capacity, and significant wave height, by direct multiplication of these parameters by the ramp function at each time step during the ramp period. Boundary conditions are also slowly transitioned from the initial condition by direct multiplication of the boundary conditions by the ramp function at each timestep during the ramp period. The ramp period is usually on the order of a few hours to a few days and should not be confused with the spinup period. The spinup period is the time it takes for the effects of the initial condition to disappear and can be much longer than the ramp period.
Hydrodynamics
Coupling of Velocity and Water Level – SIMPLEC Algorithm
The governing equations are solved in a segregated manner in which each governing equation is solved separately in a sequential manner within an iteration loop in order to obtain a converged solution. Coupling between the velocity (momentum) and water level (continuity) is achieved with the SIMPLEC algorithm (van Doormal and Raithby 1984). The main difficulty in solving the momentum equations is that the water level is not known a priori and must be calculated as part of the solution. The solution algorithm procedure is described below. First, an initial pressure is estimated based on the previous timestep value. Then, the momentum equations are solved for the corresponding velocity:
(28) 
where includes all the remaining terms. The cell face velocities are calculated with a Rhie and Chow (1983) type interpolation method
(29) 
where , is the implicit relaxation coefficient (set to 0.8 by default), and denotes linear interpolation (see Implicit Relaxation of the present Chapter). This method avoids the checkerboard oscillations associated with the collocated grid. It is noted that this approach is slightly different from that originally proposed by Rhie and Chow(1983) and is used by several authors such as Lai (2010). The current approach was found to be significantly more stable.
Next, the velocity and pressure corrections are defined such that both the momentum and continuity equations are satisfied:
(30a,b) 
In the SIMPLEC algorithm, the velocity correction is related to the pressure correction by
(31) 
where
Using and substituting in the continuity equation yields the semidiscrete water level correction equation
(32) 
Note that at convergence, and the above equation reduces to the continuity equation. Once the pressure correction equation is solved, the cellcentered water levels and current velocities are corrected. The cell face velocities are also corrected as , in which the velocity correction is given by .
Summary of the SIMPLEC Algorithm:
 Guess the water level and pressure field
 Solve the momentum equations (Equation 28) to obtain
 Use the Rhie and Chow (1983) momentum to determine the velocities and fluxes at cell faces (Equaton 29)
 Solve the pressure equation (Equation 34) to obtain
 Use the correction equations to adjust the velocities and water levels
 Treat the corrected water level and pressure field as a new guess, and repeat this procedure from Step 2 until convergence is achieved.
Wetting and Drying
During numerical simulations of the surface water flows with sloped beaches, sand bars, and islands, the landwater interface changes with time. This means that it is possible for nodes at the landwater interface to be wet or dry throughout a given simulation. In CMS, a threshold water depth is used to judge drying and wetting. If the depth at the cell center is larger than the threshold value (i.e. a small value such as 0.02 m for field cases), then the node is considered to be wet. Similarly, if the depth at the cell center is smaller than the threshold value, then the node is considered to be dry. For the implicit solver, all of the wet and dry cells are included in the matrix solver. Dry cells are assigned with a zero velocity.
Salinity Transport
Transport Equations
The sediment transport equations are discretized using the methods described in the previous section entitled General Transport Equation and are not repeated here.
Laplace Equation
The option is provided to estimate the initial salinity based on the solution of a 2DH Laplace equation with given boundary conditions and scattered salinity values. The Laplace equation is discretized similarly to the transport equation as
(33) 
where
= outward normal gradient of at cell face (f)
length of cell face (f) [m].
Sediment Transport and Morphology Change
The socalled semicoupled sediment transport model proposed by Wu (2004) is adopted here, in which the sediment calculations are decoupled from the hydrodynamic calculations, but the sediment transport, bed change, and bed material sorting equations are simultaneously solved in a coupled form at each timestep.
Transport Equations
The sediment transport equations are discretized using the methods described in the previous section entitled General Transport Equation and are not repeated here.
Bed Change Equations
The fractional bed change equation is discretized as
(34) 
where:
= morphologic acceleration factor
= bed slope coefficient []
= magnitude of the fraction bed load at the cell face f calculated from previous iteration [kg/m/s]
= bed slope calculated at the cell face.
The purpose of the morphologic acceleration factor is to speed up the bed change so that the simulation time represents approximately the change that would occur in . This factor should be used with caution and only for idealized cases or time periods which are periodic (mainly tidal). If timevarying winds or waves are important processes for driving sediment transport, then it is recommended to use reduced or representative wind and wave conditions. Since the CMS runs relatively fast, it is generally recommended not to use the morphologic acceleration factor when validating the sediment transport model using hindcast measurements. The morphologic acceleration factor is useful, however, when simulating idealized cases or analyzing project alternatives.
Bed Material Sorting
The bed material sorting equation (Equation 53) is discretized as
(35) 
where is the change in the top elevation of the second bed layer, and for and for . The mixing layer or active layer thickness calculation is slightly modified to avoid excessively small layers as
(36) 
where is the bed form height and and are the user specified the minimum and maximum mixing layer thicknesses, respectively.
The thickness of the second layer is calculated as . The bed material gradation in the second layer is calculated from the following discretized form of Equation 54:
(37) 
In order to avoid sediment layers from becoming extremely thin or thick, a layer merging and splitting algorithm is implemented between layers 2 and 3. Here, the subscript 2 corresponds to the second layer. To illustrate the bed layering process, Figure 38 shows an example of the temporal evolution of seven bed layers during erosional and depositional regimes.
It is noted that in order to maintain a constant number of bed layers, the bottom two layers are merged if the second layer is split or a new layer is created at the bottom using the bed composition and thickness of the bottom layer.
Avalanching
When the slope of a noncohesive bed, , is larger than the angle of repose, , the bed material will slide (avalanche) to form a new slope approximately equal to the angle of repose.
The process of avalanching is simulated by enforcing while maintaining mass continuity between adjacent cells. The following equation for bed change due to avalanching is obtained by combining the equation for angle of repose and the continuity equation between two adjacent cells:
(38) 
where is the cell center distance between cells P and N , is the cell area, is an underrelaxation factor (approximately 0.250.5), and H(X) is the Heaviside stepfunction representing the activation of avalanching and equal to 1 for and 0 for X < 0. The sign function, (sgn X) , is equal to 1 for and 1 for and accounts for the fact that the bed slope may have a negative or positive sign. Equation (38) is applied by sweeping through all computational cells to calculate and then modifying the bathymetry as . Because avalanching between two cells may induce additional avalanching at neighboring cells, the above sweeping process is repeated until avalanching no longer occurs. The underrelaxation factor, , is used to stabilize the avalanching process and to avoid overshooting since the equation is derived considering only two adjacent cells but is summed over all (avalanching) neighboring cells. Equation 38 above may be applied to any grid geometry type (i.e. triangles, rectangles, etc.) and for situations in which neighboring cells are joined at corners without sharing a cell face.
Hard bottom
The sediment transport and bed change equations assume a loose bottom in which the bed material is available for entrainment. However, hard bottoms may be encountered in practical engineering applications where bed materials are nonerodible, such as bare rocks, carbonate reefs, and concrete coastal structures. Hardbottom cells in CMS are handled by modifying the equilibrium concentration as in both the sedment transport and bed change equations. The bedslope term in the bed change equation is also modified so that only deposition (no erosion) may occur at hardbottom cells.
Coupling of Sediment Transport, Bed Change, and Sorting Equations
For a semicoupled sediment transport model, the sediment calculations are decoupled from the hydrodynamics but the sediment transport, bed change, and bed material sorting equations are coupled at the timestep level and thus solved simultaneously. A modified form of the iteration procedure of Wu (2004) is implemented in CMS. The equations are obtained by substituting into the bed change and sorting equations and then substituting the sorting equation into the bed change equation.
Summary of Sediment Transport SemiCoupling Procedure:=
 Calculate bed roughness’s and bed shear stresses
 Calculate the mixing layer thickness
 Calculate the potential sediment concentration capacity
 Guess the new bed composition as
 Calculate the fractional concentration capacity
 Solve transport equations for each sediment size class for
 Calculate the total and fractional bed changes
 Determine the bed sorting in the mixing layer; Repeat Step 5 and iterate until convergence.
 Update the bed elevation as
 Calculate the bed gradation in the bed layers below the mixing layer
 Calculate avalanching
 Correct the sediment concentration due to flow depth change
When using the explicit timestepping scheme, the sediment transport and morphology change are solved at a timestep which may be equal to or multiples of the hydrodynamic timestep for efficiency. However, when using the implicit timestepping scheme, the timestep is relatively big and on the order of 10 min. Because the timestep is so big, it is not necessary to use different timesteps for hydrodynamics and sediment transport. In addition, using the same timestep for hydrodynamics and sediment transport provides a better mass balance. For these reasons, the sediment transport timestep is always set to the hydrodynamic timestep in the implicit timestepping scheme.
Coupling Procedure of CMSFlow and CMSWave
CMSFlow and CMSWave can be run separately or coupled together using a process called steering (Figure 39). The variables passed from CMSWave to CMSFlow are the significant wave height, peak wave period, wave direction, wave breaking dissipation, and radiation stress gradients. CMSWave uses the updated bathymetry (if sediment transport is turned on), water levels, and current velocities from CMSFlow.
The time interval at which CMSWave is run is called the steering interval. Currently, the steering interval is constant and the input spectra in CMSWave must be at constant intervals without any gaps. A schematic showing the steering process is shown in Figure 310 in which the simulations time is plotted as a function of the computational time.
For CMS versions prior to version 4.0, the steering process is controlled in the SMS interface with communication files because the CMSWave and CMSFlow models are separate executables. For CMS v4.0 and above, both CMSFlow and CMSWave are contained within a single executable (known as the inline executable) and the steering process is controlled by an interval steering module. Two main advantages of the inline executable and inline steering module are: 1) the model runs faster because there is no need to use communication files or reinitialize the models (memory allocation, variable initialization, etc.); and 2) the inline executable makes the improvement and maintenance of the steering module easier for the developers and also makes the code portable for other operating systems. The inline steering process is summarized as follows:
 CMSWave is run for the first two timesteps and the wave information is passed to CMSFlow (Figure 310). If specified, the surface roller model is run on the wave grid and the roller contributions to the radiation stresses are added to the wave radiation stresses.
 The wave height, period, dissipation, radiation stress gradients, and wave unit vectors are interpolated spatially from the CMSWave grid to the CMSFlow grid.
 CMSFlow is run until the next steering interval and wave variables are linearly interpolated throughout time during the specified steering interval. At each flow timestep, variables such as wave length and bottom orbital velocities are updated using the current water depths and current velocities.
 Water levels, current velocities, and bed elevations are estimated for the next wave timestep and are interpolated from the CMSFlow grid to the CMSWave grid.
 CMSWave is then run again for the following timestep.
 Step 25 are repeated until the end of the simulation.
Spatial Interpolation and Extrapolation
CMS allows the user to use the same or different grids for CMSFlow and CMSWave. If the same grid is used, then no spatial interpolation is carried out. If different grids are used, then spatial interpolation is necessary in order to transfer information from one model grid to another model grid. The interpolation of wave variables from the CMSWave grid to the CMSFlow grid is done using a combination of bilinear and linear triangular interpolation methods. Bilinear interpolation is applied at nonjointed cells (i.e., cells with four neighbors) and trianular interpolation at jointed cells (cells with more than four neighbors). If the extents of the CMSWave and CMSFlow grids are different (e.g. if the CMSFlow grid is smaller), then the extrapolation of variables is necessary in order to avoid boundary problems with the models. The spatial extrapolations for different variables are given by
(39) 
(40) 
(41) 
(42) 
(43) 
(44) 
where is the mean water surface elevation, is depthaveraged current velocities, is the bed elevation, is the significant wave height, is peak wave period, and is the wave unit vector. The superscripts m and n indicate the CMSWave and CMSFlow timesteps, respectively. The subscripts wave and flow indicate variables on CMSWave and CMSFlow grids respectively. The subscripts P and N indicate variables at the extrapolated and nearest neighbor cells. The extrapolation function (f_{ext}) is given by
(45) 
Here, is the distance from cell P to N, and is the extrapolation distance. The extrapolation distance is assigned to each computational grid and can be either automatically calculated by CMS or specified by the user. The mean water surface elevation, peak wave period, and wave unit vectors are extrapolated using a nearest neighbor. This approach is more physically accurate than extrapolating only to a certain distance. For example, water levels are controlled mainly by tides along the coast, and the spatial variation is usually much smaller than the tidal range. Extrapolating bed elevations from a boundary can lead to sharp changes in bathymetry in the wave model and instability problems in both the wave and flow models. A better approach is to extrapolate the bed change. It is noted that this is only for extrapolation. For internal cells the actual bed elevations are interpolated from the flow grid to the wave grid. Finally, careful attention is needed in determining the nearest neighbor so that variables are not extrapolated over inactive portions of the grid (e.g., interpolating values in a bay using values from the open ocean).
Temporal Interpolation and Prediction
Because CMSWave requires the water surface elevation at times that are ahead of the hydrodynamic model, the water surface elevation and currents must be predicted for the CMSWave timestep. If the steering is relatively small (<30 min), then the values from the last time step may be used without significant error as
(46) 
(47) 
(48) 
where superscripts m and n indicate the CMSWave and CMSFlow time steps, respectively. The subscript flow indicates the variables on the CMSFlow grid. For many coastal engineering applications, it is desirable and common to use relatively large steering intervals of 2 to 3 hours. For large steering intervals, the change in water depth has the largest influence on the nearshore wave heights. Therefore, when using large steering intervals, it is desirable to make better predictions of water levels and not use the previous timestep water levels. In cases where the relative surface gradients at any time are much smaller than the mean tidal elevation, a better approximation of water level may be obtained by decomposing the water level into
(49) 
where is the mean water level and is a variation around the mean due to due to tidal, wave, and wind generated surface gradients. The symbol can be approximated from water level boundary conditions and is generally much larger. The water level variation may be neglected or approximated as
(50) 
For most coastal inlet applications, the above equation is a better representation of the water surface elevation and is used as the default in CMS. After spatially interpolating the wave height, period, dissipation, and forcing onto the CMSFlow grid, the variables are linearly interpolated in time. The wavelength, bottom orbital velocities, and mean wavecurrent bottom friction are then updated including currentwave interactions at each timestep.