M2D Report2
List of Figures
Figure 1. | Comparison of longshore current speed values for Visser Case 4 | 12 |
Figure 2. | Definition of continuity (conservation of volume) of sediment transport | 25 |
Figure 3. | Vertical distributions of mixing coefficient due to current and wave | 29 |
Figure 4. | Schematic hard-bottom grid and potential transport rates | 35 |
Figure 5. | Corrected hard-bottom depths and adjusted transport rates | 37 |
Figure 6. | Cell and variable definitions for CMS-M2D | 40 |
Figure 7. | Control volume definition for x-momentum equation | 41 |
Figure 8. | Control volume definition for y-momentum equation | 44 |
Figure 9. | Bottom-friction coefficient adjustment versus total water depth | 49 |
Figure 10. | Control volume definition for continuity equation | 53 |
Figure 11. | Cell numbering system | 54 |
Figure 12. | CMS-M2D grid and ADCIRC mesh for Shinnecock Inlet, NY | 80 |
Figure 13. | Time discretization of a sinusoidal water-level curve | 81 |
Figure 14. | Example CMS-M2D and STWAVE grid orientations | 83 |
Figure 15. | Example CMS-M2D and STWAVE grid spacing | 84 |
Figure 16. | Grid frame around scattered data points | 87 |
Figure 17. | Map-> 2D Grid and Refine Point dialogs | 89 |
Figure 18. | SMS-generated CMS-M2D grid | 90 |
Figure 19. | Arc options for developing coastlines | 91 |
Figure 20. | CMS-M2D Boundary Conditions dialog | 93 |
Figure 21. | Extract Boundary Conditions dialog | 96 |
Figure 22. | Extract Tidal Constituents dialogs | 96 |
Figure 23. | Assign Cell Attributes dialog | 98 |
Figure 24. | Model Checker dialog | 99 |
Figure 25. | CMS-M2D Model Control: Model Control dialog | 100 |
Figure 26. | CMS-M2D Sediment Transport Option dialogs | 102 |
Figure 27. | CMS-M2D Model Control: Model Parameters dialog | 103 |
Figure 28. | CMS-M2D Model Control: Time Control dialog | 105 |
Figure 29. | CMS-M2D Model Control: Output Control dialog | 105 |
Figure 30. | CMS-M2D Model Control: Tidal Constituents dialog | 107 |
Figure 31. | Steering Module dialog | 109 |
Figure 32. | Steering Module Restart Options dialog | 109 |
Figure 33. | Configuration for channel infilling tests | 125 |
Figure 34. | Surface elevation, current velocity, and concentration calculated by AD formulation with Lund-CIRP formula | 125 |
Figure 35. | Channel profile change by AD formulation, grain size = 0.2 mm | 126 |
Figure 36. | Channel profile change by TL formulation, grain size = 0.2 mm | 126 |
Figure 37. | Channel shape comparison between AD and TL formulations, grain size = 0.2 mm | 127 |
Figure 38. | Channel profile change by AD formulation, grain size = 0.4 mm | 128 |
Figure 39. | Channel profile change by TL formulation, grain size = 0.4 mm | 128 |
Figure 40. | Bathymetry and CMS-M2D computational grid for Ocean City Inlet | 129 |
Figure 41. | Bathymetry and STWAVE computational grid for Ocean City Inlet | 129 |
Figure 42. | Regional ADCIRC mesh and bathymetry for Maryland coast | 130 |
Figure 43. | Calculated peak water-surface elevation along Maryland coast for December 1992 storm | 131 |
Figure 44. | Calculated velocity field at Ocean City Inlet during 4.8 m wave condition, December 1992 storm | 132 |
Figure 45. | Calculated velocity fields at 10-hr intervals at Ocean City Inlet, December 1992 storm | 133 |
Figure 46. | Calculated depth at Ocean City Inlet, December 1992 storm | 135 |
Figure 47. | Calculated depth change at Ocean City Inlet, December 1992 storm | 136 |
Figure 48. | Aerial photograph of Sebastian Inlet with hard bottom shown | 137 |
Figure 49. | Bathymetry and CMS-M2D computational grid for Sebastian Inlet | 138 |
Figure 50. | Bathymetry and STWAVE computational grid for Sebastian Inlet | 139 |
Figure 51. | Detail of CMS-M2D grid resolution and bathymetry at Sebastian Inlet | 139 |
Figure 52. | Hard-bottom cells specified for Sebastian Inlet | 140 |
Figure 53. | Representative calculated ebb and flood currents at Sebastian Inlet | 141 |
Figure 54. | Calculated depth change at Sebastian Inlet after 360 hr | 142 |
List of Tables
Table 1. | Features of Calculation of Pick-up and Deposition Rates | PAGEREF TabFeaturesMethodPickupDepo 76 |
Table 2. | Boundary Condition Types | 171 |
Table 3. | Criteria for Flooding and Drying Algorithms | 185 |
Table 4. | Definitions of Parameters in Grid File | 225 |
Table 5. | Cell Types and Their Function in CMS-M2D | 228 |
Table 6. | Model Control File Specifications | 231 |
Conversion Factors, Non-SI to SI Units of Measurement
|
Non-SI units of measurement used in this report can be converted to SI units as follows:
Multiply | By | To Obtain |
acres | 4,046.873 | square meter |
cubic feet | 0.02831 | cubic meters |
cubic yards | 0.7645549 | cubic meters |
feet | 0.3048 | meters |
miles (U.S. nautical) | 1.852 | kilometers |
miles (U.S. statute) | 1.609347 | kilometers |
Preface
The tidal and nearshore circulation model CMS-M2D is being developed in the Coastal Inlets Research Program (CIRP) as part of the Coastal Modeling System (CMS). The CMS is a suite of coupled models operated in the Surface-water Modeling System, which is an interactive and comprehensive graphical user interface environment for preparing model input, running models, and viewing and analyzing results. The CIRP is administered at the U.S. Army Engineer Research and Development Center (ERDC), Coastal and Hydraulics Laboratory (CHL) under the Navigation Systems Program for Headquarters, U.S. Army Corps of Engineers (HQUSACE). Mr. Barry W. Holliday was formerly HQUSACE lead technical monitor for the CIRP. Dr. Sandra K. Knight and Mr. James E. Clausner, CHL, are the Technical Director and Associate Director, respectively, for the Navigation Systems Program. Dr. Nicholas C. Kraus, Senior Scientists Group (SSG), CHL, is the CIRP Program Manager. This report is the second in a series on CMS-M2D, documenting upgrades to the circulation model and new capabilities to calculate sediment transport and morphology change at coastal inlets and adjacent beaches and waterways.
- The mission of the CIRP is to conduct applied research to improve USACE capability to manage federally maintained inlets, which are present on all coasts of the United States, including the Atlantic Ocean, Gulf of Mexico, Pacific Ocean, Great Lakes, and U.S. territories. CIRP objectives are to advance knowledge and provide quantitative predictive tools to (a) make management of coastal inlet navigation projects, principally the design, maintenance, and operation of channels and jetties, more effective and reduce the cost of dredging, and (b) preserve the adjacent beaches in a systems approach that treats the inlet and beach together. To achieve these objectives, the CIRP is organized in work units conducting research and development in hydrodynamic, sediment transport and morphology change modeling; navigation channels and adjacent beaches; inlet structures and scour; laboratory and field investigations; and technology transfer.
- This report was prepared by Dr. Adele Militello Buttolph, Coastal Analysis LLC; Dr. Christopher W. Reed, URS Corporation; Dr. Kraus, CHL; Dr. Nobuyuki Ono, presently at ECOH Corporation, Tokyo, Japan, and formerly at Kyushu University and guest researcher, CHL; Drs. Magnus Larson, Benoit Camenen, and Hans Hanson, Lund University, Sweden; Mr. Ty Wamsley, CHL, and Dr. Alan K. Zundel, Environmental Modeling Research Laboratory, Brigham Young University. Dr. Buttolph was Principal Investigator for CMS-M2D model development. Dr. Frank S. Buonaiuto, State University of New York at Stony Brook, Marine Sciences Research Institute, Mr. Kenneth J. Connell, Coastal Engineering Branch, CHL, and Mr. Mitchell E. Brown, SSG, CHL, contributed to model development and testing. Work was performed under the general administrative supervision of Mr. Thomas W. Richardson, Director, CHL, and Dr. William D. Martin, Deputy Director, CHL.
- At the time of publication of this report, Dr. James R. Houston was Director of ERDC, and COL Richard B. Jenkins, was Commander and Executive Director.
Introduction 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).
Governing 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:
where = wind drag coefficienth=still-water depth relative to a specific vertical datum | ||
h=deviation of the water-surface elevation from the still-water level | ||
t=time | ||
qx=flow per unit width parallel to the x-axis | ||
qy=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 | ||
Dx=diffusion coefficient for the x direction | ||
Dy=diffusion coefficient for the y direction | ||
f=Coriolis parameter | ||
=bottom stress parallel to the x-axis | ||
=bottom stress parallel to the y-axis | ||
=surface stress parallel to the x-axis | ||
=surface stress parallel to the y-axis | ||
=wave stress parallel to the x-axis | ||
=
|
Component velocities are related to the flow rate per unit width as:
For the situation without waves, bottom stresses are given by:
where U = total current speed and Cb = empirical bottom-stress coefficient. The total current speed is:
The bottom-friction coefficient is calculated by:
where C = Chezy coefficient given by:
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):
where = wave angular frequency, H = wave height, and k = wave number.
Surface wind stresses are given by:
where
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:
where Sxx, Sxy, and Syy 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):
where
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:
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):
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:
where describes the lateral mixing below trough level (Smith, Larson, and Kraus 1993) and is expressed as (Kraus and Larson 1991):
where = empirical coefficient representing the lateral mixing strength, and um = amplitude of the horizontal component of the wave orbital velocity at the bottom given by:
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:
where the weighting parameter is given by:
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
5 m). A short review describing the development of the variable wind-stress coefficient implemented in CMS-M2D is presented here.
Figure SEQ Ch4Fig \* MERGEFORMAT 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):
where
The friction velocity is defined in terms of the surface wind stress
as (Hsu 1988):
where Km = 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, C10, can be expressed as (Hsu 1988):
where W10 = 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):
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 C10 and W10 defined by Equations 31 and 32 are specified as Cd 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.
Governing 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:
- Watanabe (1987) total load formulation.
- Lund-CIRP (Camenen and Larson 2005, 2006) total load formulation (combined suspended and bed load).
- 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 . 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:
where
The shear stress at incipient motion is:
(34) |
where
For the case of currents only (no waves) the bed shear stress
( SEQ Ch4Eq \* MERGEFORMAT 35) |
where is the current friction factor and in Equation 33. The current friction factor is given as (van Rijn 1998):
(36) |
in which d = total depth defined as , and ksd = Nikuradse equivalent sand grain roughness obtained from:
( SEQ Ch4Eq \* MERGEFORMAT 37) |
If waves are present, the maximum bed shear stress is given by a data-based method of (Soulsby 1997):
( SEQ Ch4Eq \* MERGEFORMAT 38) |
where = mean shear stress by waves and current over a wave cycle, = wave bed shear stress, and = angle between the waves and the current. The mean wave and current shear stress is:
( SEQ Ch4Eq \* MERGEFORMAT 39) |
The wave bed shear stress is given by:
(40) |
The wave friction factor is (Nielsen 1992):
(41) |
where R is the relative roughness defined as:
(42) |
where Aw is the semi-orbital excursion defined as:
(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 ks is assumed to consist of three components, namely grain related roughness ksd, form-drag roughness ksf, and sediment-related roughness kss (Soulsby 1997). Total roughness is obtained as the linear sum of these components:
(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):
(45) |
in which Hr = ripple height, and Lr = ripple length, which are computed from:
( SEQ Ch4Eq \* MERGEFORMAT 46) |
for a current (Soulsby 1997) and from:
( SEQ Ch4Eq \* MERGEFORMAT 47) |
for waves (van Rijn 1993), where = mobility parameter that is defined by:
(48) |
and s = specific density given by:
(49) |
The sediment-related roughness is given by Wilson (1966, 1989):
( SEQ Ch4Eq \* MERGEFORMAT 50) |
where = Shields parameter for current or waves (i = c, w, respectively, for current and waves). Expressions for the Shields parameter for currentsand for wavesare given by:
(51) | |
(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 fc and waves fw are separately computed according to Soulsby (1997) and Swart (1974), respectively as:
(53) | |
(54) |
where =Von Karman constant (taken to be 0.4). Shear stress from the current is given by:
(55) |
whereas the maximum shear stress from the waves is estimated from:
(56) |
and the mean shear stress is obtained from:
(57) |
where a sinusoidal wave is assumed. The corresponding shear velocities for current and waves, respectively, are defined as:
(58) |
for the current, and:
(59) |
for the waves (based on the maximum shear stress).
Bed load
A general transport formula for bed load under combined waves and current was developed by Camenen and Larson (2005):
(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 and are the mean and maximum Shields parameter for waves and current combined, respectively, not including the bed form roughness. The coefficient an = 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 and 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:
(61) |
in which:
(62) |
and
(63) |
where Twc = period of positive flow, Twt = period of negative flow, and uw = time-varying horizontal bottom orbital velocity ( ). Analytical solutions to these equations are available for simple expressions of uw. The friction factor for combined wave and current flow is calculated from Madsen and Grant (1976):
(64) |
where and fc and fw 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:
(65) |
and the Shields parameter applied in the initiation of motion term is:
(66) |
where and are the Shields parameters for current only and waves only, respectively, and the index m denotes the mean value. For sinusoidal waves, . 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):
(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:
(68) |
where the transport qbc 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 qs 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):
(69) |
where wf = sediment fall speed, cR = reference concentration, and = sediment diffusivity (or mixing coefficient). The transport qs 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 cR and , the total shear stress should be applied, including the bed form roughness. The reference concentration is given by:
(70) |
where the coefficient is determined by the following relationship:
(71) |
and the dimensionless grain size is defined by
. The sediment diffusivity is related to the energy dissipation according a generalization of the treatment of Battjes (1975):
(72) |
where De is a total effective dissipation given by:
(73) |
in which kc, kw, and kb are coefficients. The dissipations from bottom friction due to current Dec and from bottom friction due to waves Dew are expressed as:
(74) | |
(75) |
and the dissipation due to wave breaking Deb is directly provided from a wave transformation model (e.g., STWAVE. WABED). The coefficient kb corresponds to an efficiency factor and was found to be 0.017, whereas kc and kw 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:
(76) |
where= Schmidt number and the subscript i = c or w denotes current and waves, respectively. The Schmidt number is calculated from the empirical relationships:
(77) |
where A1 = 0.7 and A2 = 3.6 for current, and A1 = 0.09 and A2 = 1.4 for waves. In the case of coexistent waves and current, a weighted Schmidt number is applied according to:
(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 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:
(79) |
in which
The suspended sediment concentration can be solved from Equation 79 at a sediment transport time step.
Bed level change is obtained from the sediment continuity equation:
(80) |
where qbx = bed load transport rate parallel to the x-axis, qby = bed load transport rate parallel to the y-axis, and p = porosity of sediment.
In Equation 80, the terms include the bed load component, and the term P-D corresponds to the suspended load component. Figure 2 shows the definition of sediment continuity in a water column, where qsx and qsy 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 qsx and qsy may be defined as follows from Equation 79:
(81) | |
(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:
(83) | |
(84) |
where c = equilibrium concentration of suspended sediment at a given elevation, and z = vertical coordinate. Both ca and c0 are reference concentrations defined at z = a. Because the upward flux of sediment depends on the bed shear stress, ca is determined from the bed shear stress calculated from the local hydrodynamic conditions. Representation of ca 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, c0 is specified from solution of Equation 79.
Assuming that the suspended concentration is in equilibrium , then the basic equation for suspended sediment concentration can be written:
(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:
(86) |
where F(z) is a function of the vertical concentration distribution. The relationship between the reference concentration c0 and the depth-averaged concentration C is:
(87) |
Thus, c0 may be written in the following form by introducing a conversion function d:
(88) |
The bed level change due to suspended load is based upon the difference of the two types of reference concentration:
(89) |
This equation implies that erosion occurs if ca > c0, and accretion occurs if ca < c0. In the AD model, three methods of specifying ca and c0 (that is, )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 ca in the van Rijn formula is given by:
(90) |
where a = reference height and s,max = maximum bed shear stress given by Equation 38.
Features 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:
(91) |
The wave-related shear stress is given as:
(92) | |
(93) |
where = roughness height defined as:
(94) |
in which ksd and kss are calculated from Equations 37 and 50, respectively.
Bed concentration in the van Rijn model is defined at the height a as:
(95) |
where Hr = ripple height. If ripples are present, the total roughness height is modified as:
(96) |
where, ksf = 20(Hr2/Lr) is specified for the van Rijn formula, and ksf = 7.5(Hr2/Lr) 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 to determine c0 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:
(97) |
where u*c = current-related bed shear velocity expressed as:
(98) |
andis a coefficient obtained from:
(99) |
Figure3. Vertical distributions of mixing coefficient due to current and waves
The wave-related mixing coefficient is:
(100) |
where
(101) | |
(102) |
in which the parameter = height from the bed given by
, H = significant wave height, and T = significant wave period. If waves and current coexist, the combined mixing coefficient is given by:
(103) |
By applying the expression for cw, the concentration profile can be derived from Equation 85, and the conversion parametercalculated. In the AD-model based on the van Rijn equations, two methods are implemented to calculate . One is based on an exponential profile employing a depth-averaged mixing coefficient, and the other uses the original van Rijn profile, where the is obtained by numerical integration of Equation 85. Assuming an exponential profile of suspended sediment concentration,is obtained analytically as:
(104) |
where
(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 cR and ε are presented in the section describing the Lund-CIRP formula.
Horizontal diffusion coefficient
Horizontal diffusion coefficients Kx and Ky for the depth-averaged sediment concentration in the AD equation are calculated by the equation proposed by Elder (1959):
(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.
(107) |
wheredescribes the lateral mixing below trough level (Smith, Larson, and Kraus, 1993) and is expressed as (Kraus and Larson 1991):
(108) |
where = 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:
(109) |
where the weighting parameteris given by:
(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:
(111) | )
The Lund-CIRP formula for bed load is described previously in this chapter. Scaling Factors for Bed Load and Suspended Load Transport RatesAdjustments 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 EquationDepth 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 dtsed and the other over which to compute averaged transport rates and morphology change dtmorph. Directional instantaneous transport rates and are calculated at intervals of dtsed and averaged over the time interval dtmorph to obtain the time-averaged transport rates and 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 and as (Watanabe 1987):
in which = empirical slope coefficient with typical range of 5 to 30, and 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 dtmorph and for total load formulations is given by:
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, dtsed and dtmorph respectively. Also, the bed load components in Equation 80 are modified by the slope effect as: where . Finally, the sediment continuity equation is solved to compute depth change at the time interval dtmorph and is rewritten as:
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 MethodologyRepresentation 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:
The general procedure for treating hard bottom is:
The correction procedure involves four steps:
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 is greater than the hard-bottom level , the outgoing transport rates are corrected to match the criterion . Over one morphologic time step, the total amount of material , expressed as a vertical quantity that may leave the cell without violating the hard-bottom constraint is:
where qin = total influx of sediment from the cell during the morphologic time step. The correction factor K is given as:
where qout = total outflux of sediment from the cell during the morphologic time step. To ensure that 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 at all cells, and the hard-bottom level is at (units) at all cells (z is positive downwards from the original bottom level at 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. Figure SEQ Ch4Fig \* MERGEFORMAT 4. Schematic hard-bottom grid and potential transport rates
Correction loop 1:
Correction loop 2:
Correction loop 3:
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. Figure SEQ Ch4Fig \* MERGEFORMAT 5. Corrected hard-bottom depths and adjusted transport rates
Avalanching MethodologyControls 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.
Numerical Implementation of Governing Equations for Water MotionGrid definition and discretization of the governing equations are described in this chapter. The governing equations of CMS-M2D are solved numerically on a staggered rectilinear grid by a finite-volume calculation approach. All spatial derivatives are approximated with central differences with the exception of the advection terms, for which an upwind algorithm is applied. The momentum equations are solved first, followed by the continuity equation wherein updated velocity values are applied. To maximize efficiency of memory for domains with complex shoreline perimeters, the CMS-M2D grid is stored in a 1D array rather than a 2D matrix. The cell numbering convention is described in this chapter. Information is also given on the Courant stability criterion, which provides the user with understanding and an initial approximation of an appropriate model time step for the explicit solution method.
Rectilinear GridThe governing equations are solved on a discretized domain where cells are defined on a rectilinear grid (which can be regular or irregular), as shown in Figure 6. Each cell has indices i and j that correspond to its position along the x- and y-axes of the grid domain, respectively. Water-surface elevation is calculated at the cell center, whereas the x- and y-components of the velocity are calculated on the left face center and bottom face center of the cell, respectively. Values of flow rate per unit width, qx and qy, are calculated at the same locations as u and v, respectively. Note that the notation for flow rate per unit width is similar to that for sediment transport rate. For the present discussion on hydrodynamics, lower case q denotes flow rate per unit width. Because the best orientation of the grid may be such that the x- and y-axes do not correspond to geographical coordinates (N-S and E-W), cells in the grid are defined in a local (or grid) coordinate system. The local coordinate system is referenced to geographic coordinates by specification of the angle between the y-axis and true north. A positive value of this angle denotes clockwise rotation and the maximum rotation angle should not exceed 45 deg to retain an orientation as close to geographic as possible. Figure 6. Cell and variable definitions for CMS-M2D The 2D array format for referencing cell locations is retained here to facilitate presentation of the finite-volume approximations. Mapping of the cell indices from a 2D array to a 1D array is conducted primarily to reduce computer memory requirements. The 1D cell indexing method is discussed below. Momentum EquationsThe finite-volume scheme for the x-momentum equation will be given first, followed by that for the y-momentum equation. Calculation of the bottom stress coefficient, wind stress, eddy viscosity coefficient, and adjustment to the wave stress for shallow water follow the finite-volume approximation of the y-momentum equation because their descriptions are common for both x- and y-momentum equations. Solution of the momentum equations is conducted to calculate the x- and y-components of velocity u and v, respectively. The explicit approach implemented in CMS-M2D solves each momentum equation for the corresponding flow rate per unit width or at the present time step k+1 from values calculated at time step k. Velocity components and are then calculated from the values of and , respectively.
x-momentum equationThe x-momentum equation is solved explicitly by a finite-volume approximation for the control volume shown in Figure 7. The control volume is indicated by the dashed line. Surface and bottom stresses (wind stress and bottom friction) are approximated on the upper and lower surfaces, and momentum fluxes are approximated at each cell face. Figure SEQ Ch4Fig \* MERGEFORMAT 7. Control volume definition for x-momentum equation The x-momentum equation in finite volume form is given by:
where
and
The x-directed current velocity is given by:
and the current speed for the x-momentum equation is:
Wave stress for the x-momentum equation is given by:
y-momentum equationThe y-momentum equation is solved explicitly with a finite volume approximation for a control volume as shown in Figure 8. The control volume is indicated by the dashed line. Surface and bottom stresses are approximated on the upper and lower surfaces and momentum fluxes are approximated at each cell face. Figure SEQ Ch4Fig \* MERGEFORMAT 8. Control volume definition for y-momentum equation The y-momentum equation is finite volume form is given by:
where and
The y-directed current velocity is given by:
and the current speed for the y-momentum equation is:
Wave stress for the y-momentum equation is given by:
Bottom-friction coefficientThe bottom-friction coefficient is given by:
where Ci,j = Chezy coefficient calculated as:
where = hydraulic radius, which is dependent on the cell dimension , in the corresponding flow direction. The hydraulic radius for a given cell is calculated by:
where = cell-face averaged depth, Pi,j = wetted perimeter of the cell, and denotes or for application of Equation 138 to the x- or y-momentum equation, respectively. For the x-momentum equation, the cell-face averaged depth is calculated as:
and for the y-momentum equation, the cell-face averaged depth is calculated by:
The wetted perimeter is taken to be equal to if a cell has no impermeable walls. If walls are present, the wetted perimeter is calculated by:
where m = number of wall boundaries on sides parallel to the direction of flow for the momentum equation being solved. The total number of wall boundaries includes those corresponding to cell i,j and cell i-1,j if the x-momentum equation is being solved, or cell i,j-1 if the y-momentum equation is being solved. For example, in the situation of flow along a one-cell-wide channel adjacent to land and aligned with the x-axis, the wetted perimeter would be calculated as at all cells in the channel. The formulation of the bottom-friction coefficient takes into account wall friction, but does not distinguish it from friction along the channel bottom. Inclusion of waves in the bottom stress is implemented by representing the combined wave and current terms of Equations 11 and 12, respectively, as:
where
Because the wave properties are scalar values, they are located at cell centers and are face-averaged for application in the bottom stress calculation. If water in a cell approaches a thin layer, numerical instability may arise and unrealistically large values of water level and current speed calculated. In a real system, water encroaching onto a dry surface is expected to experience large frictional resistance. Stability in shallow water is enhanced within CMS-M2D by adjustment of the bottom-friction coefficient by the following formulation:
where the subscript adjusted denotes the modified bottom-friction coefficient, and a is a parameter that controls the gradient of the increase of the bottom-friction coefficient as the water depth approaches zero. The bottom-friction coefficient can vary within the range to , with the adjustment scaled by the exponential expression given within Equation 156. Figure 9 shows the depth-dependent adjustment of the bottom-friction coefficient as given by the exponential formulation with the parameter a set to 10. Significant increases in the bottom-friction coefficient occur for small water depth, particularly less than 0.1 m because of the presence of large roughness elements relative to the water depth. Thus, modification of the bottom-friction coefficient from its normal value is restricted to small water depth. For depths greater than approximately 0.2 m, modification of the bottom-friction coefficient is small. Figure SEQ Ch4Fig \* MERGEFORMAT 9. Bottom-friction coefficient adjustment versus total water depth Wind stressWind-stress terms are calculated by:
where Cd is calculated by Equation 31, and is calculated from W. The density ratio is specified as 0.0012. Height of the anemometer is specified in an input file. Wind forcing in CMS-M2D can vary over time, but is homogenous over space because of anticipated local- or project-level application. Eddy viscosity coefficientTwo formulations for calculation of the eddy viscosity coefficient are implemented in CMS-M2D, one for oceanic mixing and one for surf zone mixing. If waves are not present, the oceanic coefficient of eddy viscosity is calculated as:
where
These formulations are calculated for Manning’s n > 0. For frictionless cells, specified by , with . In the presence of waves and wave breaking, the eddy viscosity coefficient is calculated as:
Wavelength l is calculated by the Eckart (1952) approximation:
where = deepwater wavelength given by:
and w is the wave frequency given by:
Wave parameters are calculated at locations corresponding to indices on Dw. If waves are present, but not breaking, then the mixing coefficient is calculated as a weighted function of the oceanic and wave mixing values as:
where the weighting parameter is:
Equations 173 and 174 are calculated and applied at locations corresponding to indices on Do and Dw. Wave stress scaling for shallow waterCMS-M2D accepts wave stresses calculated by a wave model, and those stresses can vary spatially and temporally over the computational domain. Wave stresses must be mapped from the wave model onto the CMS-M2D domain, and those stresses are computed for water depths provided to the wave model. Wave properties are typically updated over time intervals of hours, during which water levels can rise and fall owing to the tide, storm surge, and wind forcing. During a CMS-M2D simulation, the total water depth will vary as a function of the total forcing. If the water level decreases, the wave stresses can contribute disproportionately to the velocity, particularly in shallow water because the wave stresses are not updated by the wave model at every time step. For example, the water depth at a specific location may be 0.5 m for the wave calculations and at the beginning of a CMS-M2D simulation. During the CMS-M2D simulation, the water level lowers to 0.05 m, but the wave stress value was calculated for 0.5 m and may be significantly larger than it would be for a depth of 0.05 m. Thus, an ad hoc scaling factor has been implemented to adjust the wave stresses in shallow water. This scaling factor reduces the wave stress in shallow water depths to maintain realistic velocities and to promote stability. In water depths of 0.35 m or less, the wave stresses are reduced by:
where = adjusted wave stress and is applied to both ]] and if the depth criterion is met. Equation 175 scales the wave stress toward zero as the water depth decreases. Wave models typically perform calculations in a local right-handed coordinate system in which positive x is directed onshore and positive y is directed alongshore. Thus, the orientation of the local coordinate system may not be aligned with global geographic coordinates. CMS-M2D performs calculations on a local coordinate system that is referenced to global geographic coordinates by the deviation of the y-axis (given in degrees) from true north. If wave stresses are included in CMS-M2D calculations, mapping of those stresses from the wave model to the circulation model must take into consideration the different coordinate systems and grid orientations. Wave stresses calculated by STWAVE or WABED can be automatically mapped to CMS-M2D within the SMS (see Chapter 8). Continuity EquationThe continuity equation is solved explicitly for water-surface elevation by a finite-volume approximation for a control volume depicted in Figure 10. Mass fluxes are approximated at each cell face. In the staggered time differencing scheme, the most recently calculated values for the u and v velocity components are entered into the spatial derivatives. The finite-volume approximation of Equation 1 is given by:
where and all other variables have been previously defined. Figure SEQ Ch4Fig \* MERGEFORMAT 10. Control volume definition for continuity equation Cell NumberingCMS-M2D stores its numerical grid in an array rather than a matrix, which saves memory for large or irregular grid domains. The numbering scheme for the computational grid is designed to be efficient such that cells that reside in computationally inactive regions do not have to be stored. For model domains that have complex shorelines, this feature can save substantial computer memory. On the other hand, for simple grids, more memory may be required because pointing arrays are necessary for cell indexing.
The cell numbering system internal to CMS-M2D does not require row and column indexing, but instead assigns each cell a unique identification number. Flather and Heaps (1975) implemented a similar cell numbering system so that only cells included in hydrodynamic calculations were stored. Non-computational regions of the grid, such as the mainland and islands, do not require assigned cell numbers, so they do not occupy space in memory. Figure 11 illustrates the cell numbering system for a small domain that includes a non-computational region indicated by hatching. Cell identification numbers of neighbors adjacent to a particular cell are stored in a pointing matrix within the model. For example, cell 9 has neighbors, with cell identification numbers 12, 10, 4, and 8 moving from top, to right, bottom, and left, respectively. This clockwise order of referencing neighbor cells is taken as convention for descriptions in this document. Figure SEQ Ch4Fig \* MERGEFORMAT 11. Cell numbering system Locations where no cells exist, such as beyond the computational domain or in non-computational regions such as islands within the domain, are not stored in the grid as cells, but can be referenced as having a cell identification number of 0. For example, cell 11 in Figure 11 lies on the boundary of the grid domain and also resides next to a non-computational area. Neighbor-cell specifications for cell 11 would be stored as 14, 0, 6, and 0. The grid does not contain any cells with a cell number of zero; this value is used in neighbor specification to indicate that a computational cell does not exist at a particular location. Courant ConditionFor an explicit solution method, an initial estimate of the maximum time step for a grid can be calculated from the Courant number x, given by (Richtmyer and Morton 1967):
The theoretical maximum value of the Courant number is unity for stability of a linear, finite-volume hydrodynamic model, in which for shallow-water equations refers to the speed of a long wave, such as the tide. For a fixed and u, the Courant number limits the value of the hydrodynamic time step . Because CMS-M2D solves the nonlinear shallow-water equations, stability will be maintained for a Courant number less than unity. Also, for a typical coastal inlet situation, multiple forcing is possible from a superposition of sources other than or in addition to the tide. Contributions to the current can be produced by wind, surface waves, and tributary discharges. Assigning a velocity component to each of these forcings, denoted with a subscript, the Courant number is more accurately defined as:
In strong or focused-flow regions of a model domain, such as at ebb shoals and in the surf zone, the wave-driven current can be significantly stronger than the tidal current. Moreover, these are typically areas of greater spatial grid resolution. The combination of strong currents and small grid cells limits the allowable size of the time step in an explicit solution scheme. Numerical Implementation of Governing Equations for Sediment Motion
Three choices are available for computing sediment transport with CMS-M2D: the Watanabe (1987) formulation, the Lund-CIRP formula (Camenen and Larson 2005, 2006), and the AD equation. The Watanabe transport formula produces a total transport rate, estimating the contributions of combined suspended sediment load and total load in a "calculation-at-a-point" formulation. The Lund-CIRP transport formula computes both bed load and suspended load transport rates, which are combined to obtain a total load transport rate. The AD equation is solved for both suspended load and bed load transport rates.
Morphology change is computed by means of two time steps, a transport rate time step and a morphology change time step. Instantaneous transport rates are computed at the transport rate time step and averaged over the morphology change time step. Averaged transport rates are then applied in the sediment continuity equation. Transport rates calculated by the Watanabe, Lund-CIRP total load formula, and AD equation are instantaneous values.
Watanabe FormulationImplementation of the Watanabe formulation computes transport rates on cell faces with the x-directed rate located on the left face center and the y-directed rate located on the bottom face center. Thus, for a given cell, and are calculated at the same locations as u and v, respectively (Figure 6). Time-averaged face-centered transport rates are applied in the sediment continuity equation to calculated depth change. Calculation of x- and y-directed transport rates and , respectively, by the Watanabe formulation is conducted in CMS-M2D as:
The bottom stresses in the situation in which no waves are present are calculated in the x and y directions as:
and
in which the velocities and are calculated as:
and
The current friction factors for the x and y directions, respectively, are:
and
If waves are present, the bottom shear stresses in the x and y directions are calculated as:
and
Bottom shear stresses calculated by Equations 187 and 188 are entered as values of and in Equations 179 and 180, respectively. Combined wave and current shear stresses for the x and y directions are given by:
and
The wave bed shear stresses for the x and y directions, respectively, are given by:
and
where the wave friction factor for the x direction is given as:
and for the y direction is:
Relative roughness for the x and y directions are given by:
and
Wave orbital velocity amplitudes for the x and y directions, respectively, for non-breaking waves are:
and
For breaking waves, the wave orbital velocity amplitudes for the x and y directions, respectively, are:
and
Lund-CIRP FormulationThe Lund-CIRP formulation computes sediment transport properties and transport rates at the cell center, which corresponds to the same location in the cell as the water-surface elevation is calculated (Figure 6). After transport rates are calculated by the Lund-CIRP formulas for all cells in the domain, they are mapped back to the cell faces such that the x-directed rate is located on the left face center and the y-directed rate is located on the bottom face center. Time-averaged face-centered transport rates are applied in the sediment continuity equation to calculated depth change. AD EquationAD equation for sediment concentrationThe AD equation for sediment concentration is solved explicitly for depth-averaged concentration by a finite-volume approximation for a control volume. Depth-averaged concentration C, pick-up rate P, and deposition rate D are approximated at each cell center, and mass fluxes are approximated at each cell face. The finite-volume approximation of Equation 79 is given by:
where
and
Pick-up rate and deposition ratePick-up rate and deposition rate can be calculated in CMS-M2D by either the van Rijn or Lund-CIRP formulation. The calculation method of P and D for van Rijn formula is described first, then that for Lund-CIRP formula. Van Rijn model (exponential profile or van Rijn profile) The pick-up rate P and for van Rijn model is given by:
Maximum bottom shear stress s,max is given by:
where cs, ws, and ms are shear stresses owing to current, wave, and combined wave and current, and are given by:
The current velocity Uc and wave orbital velocity Uw are defined at each cell center.
and wave friction factor is:
where
Calculations with CMS-M2D have shown that the van Rijn (1985) reference concentration tends to give an overestimation under strong shear stress, such as in the sheet flow region. Therefore, the following limitation is applied to reduce the overestimate:
where R is the reduction coefficient expressed as:
in which cb is the extrapolated bed concentration calculated by Equation 207 with a = d50, which reduces to:
The deposition rate D for the van Rijn model is given by:
where d is the conversion parameter that depends on the user-selected concentration profile. If an exponential profile is selected, d is calculated by using depth-averaged mixing coefficient and given by:
and the reference height a is given by:
where Hr is the ripple height. If the van Rijn profile is selected, d has to be calculated numerically because the mixing coefficient varies with height from the bed. In the AD-model with the van Rijn profile, the mixing coefficients are calculated at 20 points in the vertical direction. Mixing coefficients for the current c and for waves w are calculated, and then combined as follows:
where
With the mixing coefficient cw(z) available, the vertical distribution of the concentration c(z) can be solved by the fourth-order Runge-Kutta method. Then, d values for the van Rijn profile are calculated by numerical integration of the concentration profile:
Lund-CIRP model The pick-up rate P for the Lund-CIRP model is given by:
where ca is the reference concentration given by:
where AcR = 3.5 10-3exp(-0.3d*), and b = 4.5. Bed shear stresses due to current, wave, and combined current and wave are given by,
Friction factors for the current and for waves are given, respectively, by:
The deposition rate D for the Lund-CIRP model is given by:
The concentration profile for the Lund-CIRP formula is approximated by an exponential profile. The conversion parameter d is given by:
where the mixing coefficient is calculated as:
De is the energy dissipation rate owing to current, wave, and wave breaking:
where
and Db is given from wave model output, and and Sediment diffusion coefficientDiffusion coefficients Kx and Ky are defined at cell faces and calculated as:
Bed load transport rateBed load transport rates qbx and qby are calculated at each cell center. For the van Rijn model, they are given by:
where
Bed load transport rates for the Lund-CIRP model are given by:
where an = 12 and b = 4.5. Boundary conditionsBoundary conditions for the AD equation are given for wall boundary and for open boundary, respectively. If a cell face is defined as a wall, the suspended load flux on the cell face is set to zero.
For the cell on the open/forcing boundary, the concentration is set as the equilibrium concentration,
where is the depth-averaged concentration on the open boundary cell. Also, the gradient of suspended load flux is set as zero ( or ) for the open boundary.
Criterion of time step for AD modelCalculation of the AD equation is computed at a time step of tsed. To maintain stable calculations in an explicit solution scheme, the following condition must be satisfied:
where s is the cell size, and u is the local current speed. During calculation of the AD equation, CMS-M2D evaluates whether this criterion is met. If the criterion is not met, CMS-M2D will reduce tsed to maintain stability. Hard-bottom methodology for AD modelModification of sediment transport rate owing to hard bottom is also computed in the AD model. Hard bottom is treated separately for suspended and bed load. For suspended load, the methodology assumes that upward sediment flux will be zero if the hard bottom is exposed. The criterion for the maximum pick-up rate is expressed as:
where hhb is the hard-bottom level, and p is the porosity. For bed load, the methodology of treating hard bottom as described in Chapter 3 is applied. Sediment Continuity EquationCalculation of depth change is computed by solution of the sediment continuity equation. Before solving the sediment continuity equation, time-averaged sediment transport rates are computed from the instantaneous transport rate values. Transport rates applied in the sediment continuity calculations are total load transport rates when total load formulations are applied, or bed load transport rates when the AD equation is applied. Thus, for the sediment continuity equation, the subscript t denotes total load for total load formulations and bed load for the AD equation. Time-averaged transport rates for the x and y directions, respectively, are computed by:
in which Nm = number of instantaneous transport rate values, and l = index representing the time at which the instantaneous rates were calculated. Instantaneous rates are calculated at the sediment transport time step , and time-averaged rates apply instantaneous rates averaged over the morphologic change time step . Adjustment of the time-averaged transport rate for bottom slope is computed for the x and y directions as:
where the total time-averaged transport rates and are calculated as:
and
Depth change over the time interval is calculated as:
Updated depth is then computed by:
After new depth values are calculated for all cells in the computational domain, they are provided to the hydrodynamic portion of the model so that it can respond to changes in morphology. For the AD model, calculation of depth change is computed by a method similar to that described above, except that suspended load and bed load components are treated separately. Before solving the sediment continuity equation, time-averaged sediment transport rates are computed from the instantaneous transport rate values. Time-averaged pick-up rate and deposition rate for the suspended load component are computed by:
Time-averaged bed load transport rates for the x and y directions, respectively, are computed by:
Bed load transport rates for the AD model are defined at each cell center. Solution of the sediment continuity equation requires that the bed load transport rates are defined at the centers of cell faces. Therefore, bed load transport rates are computed on the cell faces by:
Further, adjustment of the time-averaged bed load rate for bottom slope is computed for the x and y directions as:
Depth change over the time interval is calculated as:
Updated depth is then computed by Equation 271. Boundary ConditionsSix types of boundary conditions are implemented in CMS-M2D and are distinguished as specifying a forcing or a non-forcing boundary. Table 2 lists the boundary condition types and contains a short description of each. Each boundary condition type is described next.
Water-Level Forcing Boundary ConditionsWater-level forcing can be specified through either tidal constituents or by time series of water level read from a file. These two boundary conditions apply the same calculation approach, with the distinguishing feature between them being the source of the water level prescribed at the boundary. If a cell is designated as a water-surface elevation boundary cell, the continuity equation for that cell is not solved, and the water level is updated each time step according to:
where is a prescribed value. If tidal forcing is specified, then
Tidal-constituent forcingTidal-constituent forcing in CMS-M2D can accommodate as many as eight components of the astronomical tide. The user can select from the following constituents: M2, S2, N2, K2, K1, O1, M4, and M6. Although the M4 and M6 constituents are harmonics, they can be generated in continental shelf areas and may be present as components of the tidal forcing at the grid boundary. Specification of tidal constituents requires input information for the local amplitude and phase values. Tidal-constituent forcing by CMS-M2D is simple and does not include advanced tidal calculations or relations to Greenwich tidal parameters (Schureman 1924). Tidal elevation is given by:
where A = amplitude (m), s = speed (deg/hour), and f = phase (deg) of the tidal constituent. Constituents can be eliminated from the tidal forcing by setting their amplitudes to zero. If the tidal elevation boundary condition is invoked, water levels at tidal boundary cells are computed at every time step. Tidal phase and amplitude is held constant across the tidal-constituent boundary. For applications requiring amplitude and phase variation over an open boundary, two approaches to forcing the boundary are available. Time series of tidal curves along the boundary can be developed within the SMS in which the tidal information is obtained from a tidal-constituent database. Alternatively, solutions from regional models can be applied to map water-surface elevations as boundary conditions for CMS-M2D.
Time series water-level forcingThe time series water-level forcing boundary condition is specified by water levels ( ) read into the model from an input file. Time increments of the forcing data are not required to be equal to, or in multiple increments of, the model time step. Because the model will usually have a hydrodynamic time step smaller than the time interval between water-level forcing input, CMS-M2D must be instructed on treatment of the water-level boundary value between times that its value is specified in the input file. One option is to linearly interpolate in time between input water-level values, and this is the preferred option. A second option is to keep the water level constant between input water-level values. The boundary condition input file for water level contains a parameter that sets the interpolation option.
Time Series Water-Level and Velocity Forcing Boundary ConditionThe time series water-level and velocity forcing boundary condition applies velocities at cell boundaries in addition to water levels. Application of water levels for this boundary type is identical to that previously described for water-level forcing boundary conditions. Velocities prescribed on the boundaries enter the momentum equation through the advective terms. Discretization of the momentum equations for applying velocities on the boundaries is identical to the implementation for non-boundary cells described in Chapter 4.
Time Series Flow Rate Forcing Boundary ConditionThe time series flow rate forcing boundary condition can be assigned at locations where the flow rate is known. Typical examples where this boundary condition is applied are natural tributaries (rivers, streams, and creeks), controlled-discharge canals, and anthropogenic intake or discharge locations such as for power plant cooling water or treatment plant outfalls. A flow rate boundary condition specifies discharge at one or more cells on the grid boundary. If a cell is designated as a flow rate boundary cell, the momentum equation for the appropriate cell face is not solved, and the flow-component normal to the face is prescribed. Because CMS-M2D operates on a staggered grid, virtual cells are added to the right or top of the grid at flow rate forcing boundary cells prescribed on those sides. Addition of virtual cells, where water-surface elevation and current are not calculated, is required to calculate fluxes at cell faces on the right or top of the grid. For a flow rate forcing cell located on the right side of the grid at location i+1,j, the flux is given by:
where = input flow rate directed parallel to the x-axis, and denotes an imaginary cell adjacent and to the right of the ith cell. For a flow rate forcing cell located on the left side of the grid at location
For a flow rate forcing cell located on the top side of the grid at location i,j+1, the flux is given by:
where = input flow rate directed parallel to the y-axis. For a flow rate forcing cell located on the bottom side of the grid at location
Closed, Reflective Boundary ConditionThe closed, reflective boundary condition is specified at impermeable boundaries, such as at the land-water interface. These boundaries behave as walls; water can flow parallel to the wall face, but not through it. The normal velocity at a closed boundary on the left cell face is:
and the normal velocity at a closed boundary on the bottom face is:
Closed boundaries on the right and top cell faces are prescribed at the imaginary cells i+1,j and i,j+1, respectively.
Wave-Adjusted Water-Level and Velocity ForcingSimulations for nearshore areas may have available boundary forcing information that does not include the presence of waves. The wave-adjusted water-level and velocity boundary condition modifies the prescribed forcing values to account for waves. If water-level forcing or mixed water-level and velocity forcing are applied together with radiation-stress forcing, then in this boundary condition, boundary values are adjusted for wave-driven current and set up or set down. Boundary adjustment is conducted by creating a surface deflection (set up and set down) consistent with that of the interior solution so that radiation stress-forced flow can move smoothly through the boundary (Reed and Militello 2005). It is believed that this boundary condition is new to coastal circulation modeling. Implementation of the wave-adjusted boundary conditions within CMS-M2D is automated. Thus, if radiation-stress gradients are applied as forcing, CMS-M2D will automatically apply the wave-adjusted boundary condition to adjust the prescribed forcing for the presence of waves. Adjustment to the prescribed boundary conditions is conducted by solving the 1D continuity and momentum equations normal to the boundary, in which the radiation stress gradient is included in the momentum equation. The forms of the 1D continuity and momentum equations for the x-direction are:
in which:
where h is the applied boundary condition, is the prescribed boundary value (read from a file), and are the solutions from Equations 290 and 291, and is given by:
where and denote face-centered averages that place h and at the same position as . If the boundary is specified a mixed water-surface elevation and velocity boundary, then the u-component of velocity assigned at the boundary is modified by:
where u is the applied velocity at the boundary and ub is the prescribed velocity (read from a file). The 1D continuity and momentum equations for the y-direction are:
where and are the solutions from Equations 295 and 296, and is given by:
where and are face-centered averages that place h and at the same position as . If the boundary is specified as a mixed water-surface elevation and velocity boundary, then the v component of velocity assigned at the boundary is modified by:
where v is the applied velocity at the boundary and vb is the prescribed velocity (read from a file). The wave-adjusted boundary condition is demonstrated in Appendix G for an idealized beach.
Model FeaturesCapabilities have been implemented into CMS-M2D to represent flooding and drying, and other features designed to promote flexibility and convenience of the modeling system. This chapter describes these features and provides information on model options. Topics covered are: flooding and drying, hot-start options, model spin-up, coupling with larger scale circulation models, coupling with wave models, and output options.
Flooding and DryingNumerical modeling of shallow-water areas such as beaches, flood shoals, or expansive tidal or wind-tidal flats is most accurately performed if realistic and robust flooding and drying algorithms are invoked. Representation of flooding and drying also eliminates the need to modify grid depth to avoid artificially induced instabilities. Numerical simulation of flooding and drying of cells or elements in hydrodynamic models has been performed with various techniques such as application of barrier and weir equations (Reid and Bodine 1968; Butler 1978; Wanstrath 1978) and moving boundaries (Yeh and Chou 1978). A simple and robust flooding and drying algorithm has been implemented to represent inundation and drying of shallow areas of the model domain. Criteria were established for the flooding and drying algorithms implemented in CMS-M2D and are given in Table 3. The criteria were specified so that the algorithms would produce realistic behavior, and numerical problems, such as oscillations and flutter, would be minimized.
Every wet cell is checked to determine whether or not it has dried after new water-level and velocity values have been calculated over the model domain. The criterion for a cell to be dry is:
where = total water depth, and = depth below which drying is assumed to occur. If a cell or its adjacent neighbor to the south or west has become dry, restrictions on flow between the two cells are imposed. Restrictions governing flow across the two faces are imposed for the following conditions:
and
This restriction sets the velocity at the face to zero if the velocity at the previous time step indicates transfer of water from a neighboring dry cell to the wet cell i,j. Water is allowed to flow from the wet cell to the neighboring dry cell.
This restriction sets the velocity at the face to zero if the velocity at the previous time step indicates transfer of water from dry cell i,j a neighboring wet cell. Water is allowed to flow from the neighboring wet cell to the dry cell i,j. Hot-Start CapabilitiesCMS-M2D has two hot-start capabilities, and both can be invoked during a simulation. The two hot-start types are:
An example sedparams.dat file for the Lund-CIRP total load formulation option is provided in Figure A12. The first line contains a value of 2, which informs CMS-M2D that the Lund-CIRP total load formulation is to be applied. The remaining lines contain the input options and parameters required by the Lund-CIRP formulation and morphology change calculations.
An example sedparams.dat file for the AD equation option is provided in Figure A13. The first line contains a value of 3, which informs CMS-M2D that the AD equation is to be applied. The remaining lines contain the input options and parameters required by the AD equation and morphology change calculations. In this example, a value of 3 for the suspended concentration profile denotes that the Lund-CIRP formulation will be applied. A value of 1 for the suspended concentration profile denotes that the exponential model will be applied and a value of 2 denotes that the van Rijn model will be applied.
Appendix B |
This appendix contains a listing of publications known to the authors to date that document M2D and its applications. Citations are listed in reverse chronological order.
Buttolph, A. M., Reed, C. W., Kraus, N. C., Ono, N., Larson, M., Camenen, B., Hanson, H., Wamsley, T., and Zundel, A. K. (2006). "Two-dimensional depth-averaged circulation model CMS-M2D: Version 3.0, Report 2: Sediment Transport and Morphology Change," ERDC/CHL TR-06-__, U.S. Army Engineer Research and Development Center, Vicksburg, MS (this volume).
Lin, L., Mase, H., Yamada, F., and Demirbilek, Z. (2006). "Wave-Action Balance Equation Diffraction (WABED) Model; Part 1: Tests of wave diffraction and reflection at inlets," Coastal Inlet Research Program Technical Note ERDC/CHL CHETN-III-73, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Offshore & Coastal Technologies, Inc. (2005). "An assessment of jetty shortening alternatives for Goldsmith Inlet, bay and adjacent shorelines," Project report for Town of Southold, NY, 66 pp.
Offshore & Coastal Technologies, Inc. (2005). "Ocean City hot spot study: Reassessment of design criteria and approaches," Engineering Appendix for U.S. Army Engineer District, Baltimore, General Re-evaluation Report, 109 pp.
Reed, C. W., and Militello, A. (2005). "Wave-adjusted boundary condition for longshore current in finite-volume circulation models," Ocean Engineering 32, 2121-2134
Hanson, H., and Militello, A. (2005). "Representation of non-erodible (hard) bottom in two-dimensional morphology change models," Coastal and Hydraulics Engineering Technical Note ERDC/CHL CHETN-IV-63, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Bounaiuto, F., and Militello, A. (2004). "Coupled circulation, wave, and sediment transport modeling for calculation of morphology change, Shinnecock Inlet, New York," Proceedings 8th 'International Estuarine and Coastal Modeling Conference, ASCE, 819-838.
Lin, L., Kraus, N. C, and Barcak, R. G. (2004). "Modeling sediment transport at the Mouth of the Colorado River, Texas," Proceedings 8th 'International Estuarine and Coastal Modeling Conference, ASCE, 988-1006.
Militello, A., Reed, C. W., Zundel, A. K., and Kraus, N. C. (2004). "Two-dimensional depth-averaged circulation model M2D: Version 2.0, Report 1: Documentation and user’s guide," ERDC/CHL TR-04-02, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Militello, A., and Zundel, A. K. (2003) "SMS Steering Module for coupling waves and currents, 2. M2D and STWAVE," Coastal and Hydraulic Engineering Technical Note ERDC/CHL CHETN-IV-60, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Militello, A., and Kraus, N. C. (2003). "Numerical simulation of sediment pathways at an idealized inlet and ebb shoal," Coastal Sediments 2003, CD ROM: Numerical simulation of sediment pathways at an id.pdf.
Militello, A., Parry, R. M., and Arden, H. T. (2003). "Numerical simulation of navigation channel infilling, Bay Center, Washington," Proceedings Dredging 2002 Conference, ASCE, CD ROM 40680-034-003, ISBN 0-7844-0680-4.
Militello, A., and Zundel, A. K. (2002). "Coupling of regional and local circulation models ADCIRC and M2D," Coastal and Hydraulics Engineering Technical Note ERDC/CHL CHETN-IV-42, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Militello, A., and Zundel, A. K. (2002). "Automated coupling of regional and local circulation models through boundary condition specification," Proceedings 5th 'International Conference on Hydroinformatics, International Association of Hydraulic Research, IWA Publishing, London, 156-161.
Militello, A. (2002). "Willapa entrance and Bay Center modeling," in "Study of navigation channel feasibility, Willapa Bay, Report 2," N. C. Kraus, H. T. Arden, and D. P. Simpson, eds. Technical Report ERDC-CHL-TR-00-6, U.S. Army Engineer Research and Development Center, Vicksburg, MS.
Militello, A., and Kraus, N. C. (2001). "Generation of harmonics by sea breeze in nontidal water bodies," Journal of Physical Oceanography 31(6), 1,639-1,647.
Militello, A. (2000). "Hydrodynamic modeling of a sea-breeze dominated shallow embayment, Baffin Bay, Texas," Proceedings 6th 'International Estuarine and Coastal Modeling Conference, ASCE, 795-810.
Kraus, N. C., and Militello, A. (1999). "Hydraulic study of multiple inlet system: East Matagorda Bay, Texas," Journal of Hydraulic Engineering 125(3), 224-232.
Militello, A., and Kraus, N. C. (1998). "Numerical simulation of hydrodynamics for proposed inlet, East Matagorda Bay, Texas," Proceedings 5th 'International Estuarine and Coastal Modeling Conference, ASCE, 805-818.
Militello, A. (1998). Hydrodynamics of wind-dominated, shallow embayments. Ph.D. dissertation, Division of Marine and Enviromental Systems, Florida Institute of Technology, Melbourne, FL.
Brown, C. A., and Militello, A. (1997). "Packery Channel feasibility study: Circulation and water level; Report 2 of a two-part series," TAMU-CC-CBI-96-07, Conrad Blucher Institute for Surveying and Science, Texas A&M University-Corpus Christi, Corpus Christi, TX.
Kraus, N. C., and Militello, A. (1996). "Hydraulic feasibility of the proposed Southwest Cut, East Matagorda Bay, Texas," TAMU-CC-CBI-96-03, Conrad Blucher Institute for Surveying and Science, Texas A&M University-Corpus Christi, Corpus Christi, TX.
Militello, A., and Kraus, N. C. (1995). "Investigation of the current and sediment movement in the Lower Laguna Madre, Texas, and implications for dredging in the Gulf Intracoastal Waterway," Proceedings Western Dredging Association 16th 'Technical Conference, Center for Dredging Studies, Texas A&M University, College Station, TX, 143-162.
Brown, C. A., Kraus, N. C., and Militello A. (1995). "Hydrodynamic modeling for assessing engineering alternatives for elevating the Kennedy Causeway, Corpus Christi, Texas," Proceedings 4th 'International Estuarine and Coastal Modeling Conference, ASCE, 681-694.
Brown, C. A., Militello, A., and Kraus, N. C. (1995). "Hydrodynamic assessment for elevation of the JFK Causeway, Corpus Christi, Texas," Proceedings Texas Water ’95, Texas Section of the American Society of Civil Engineers, ASCE Press, NY, 31-41.
Brown, C. A., Militello, A., and Kraus, N. C. (1995). "Hydrodynamic assessment for elevation of the JFK Causeway, Corpus Christi, Texas," TAMU-CC-CBI-95-07, Conrad Blucher Institute for Surveying and Science, Texas A&M Univeristy-Corpus Christi, Corpus Christi, TX.
Militello, A., and Kraus, N. C. (1995). "Numerical simulation of the circulation in the Lower Laguna Madre, Texas, and implications for dredging of the Gulf Intracoastal Waterway," Estuarine Research Federation 13th 'Biennial International Conference, Abstracts, 88.
Militello, A., and Kraus, N. C. (1994). "Reconnaissance investigation of the current and sediment movement in the Lower Laguna Madre between Port Isabel and Port Mansfield, Texas," TAMU-CC-CBI-94-04, Conrad Blucher Institute for Surveying and Science, Texas A&M University-Corpus Christi, Corpus Christi, TX.
Appendix C
Slosh Tests
This appendix describes the performance of CMS-M2D for slosh tests. This set of numerical integrity tests evaluates mass conservation and numerical damping properties of the model. Approximations to the full equations of motion truncate terms to some order in a manner that depends on the numerical solution method. Truncation of terms is known to introduce numerical viscosity (Roache 1976),[1] which is manifested as damping of a signal. The slosh test provides a means of examining the damping properties of numerical computation methods. In a slosh test, the initial condition of a tilted water surface is imposed, and the basin is allowed to freely propagate the resulting wave through reflection at both ends. External forcing and bottom frictional stress are set to zero so that there is no masking of the free-surface oscillations. Time-series of water elevation over the duration of the simulation exhibit model damping properties. Tests described here were patterned after those by Cialone (1989).
The slosh test was performed to examine damping properties and mass conservation for two grid configurations. Both grid configurations were rectangular, with the grid being longer parallel to the x-axis than parallel to the y-axis. A test consisted of applying a tilted water surface as the initial condition, where the direction of tilt was parallel to the x-axis. The initial deviation from still-water depth was symmetric about the channel center and increased or decreased sinusoidally along the main channel axis (Figure C1).
No variation in water level was specified across the width (parallel to the y-axis) of the channel. Initial current velocity was set to zero, and no boundary forcing, bottom friction, or wind stress was applied. The Coriolis term and lateral mixing were set to zero. Simulation duration was set to 30 days (720 hr) and the time step was set to 5 s. General grid configurations and inclusion of advective terms are given in Table C1. Conducting identical tests with and without inclusion of the advective terms allows isolation of their contribution to damping and mass conservation.
Figure C SEQ AppAFig \* MERGEFORMAT 1. Slosh test initial condition
Table C Slosh Test Configurations | ||
Test | Bottom | Advection |
1 | Flat | No |
2 | Flat | Yes |
3 | Sloped | No |
4 | Sloped | Yes |
Mass conservation was calculated as percent change in water volume over the numerical grid between time t = 0 and the end of the simulation period. Percent change in water volume
was calculated as:
(C SEQ AppAEq \* MERGEFORMAT 1) |
where
is the final volume,
is the initial volume, the subscript i indicates cell number, and n is the total number of cells in the grid.
Tests 1 and 2 were conducted to evaluate damping and mass conservation of the model for a flat-bottomed, constant-spaced grid, with and without the advective terms included in the calculations, respectively. Grid configurations and numerical parameters for Tests 1 and 2 are given in Table C2, and the grid is shown in Figure C2. Tests 1 and 2 resulted in 5 ´ 10-6 and 1 ´ 10-5 percent changes in volume, respectively.
Table C2 Slosh Test Grid Configurations and Numerical Parameters for Tests 1 and 2 | ||
Numerical Parameter | Value | |
D'''s''', m | 500 | |
D'''t''', s | 5 | |
h, 'm' | 10 | |
Length, m | 10,000 | |
Width, m | 2,500 |
Figure C2. Computational grid for Tests 1 and 2
Damping properties for Tests 1 and 2 are shown by time series of water level in Figures C3 and C4, respectively. Water levels shown for each plot were taken 0.5 Dx from the boundary and centered with respect to channel width. Time series shown for all other tests are also taken at the same cell. Axis breaks (x‑axis) were inserted into the plots so that time series at the beginning and ending portions of the simulations could be easily viewed. Test 1 does not exhibit damping over the duration of the simulation. Weak nonlinear motion owes to the nonlinear form of the continuity equation. Test 2 exhibits weak damping over the 1-month simulation. Development of nonlinearities is slightly stronger than for Test 1 owing to the contribution by the advective terms. Tests 1 and 2 were repeated with an initial surface deviation 10 times smaller than described here (results not shown). Motion induced by this smaller water-surface deviation was not sufficient to develop nonlinearities and the resulting slosh motion for both cases (with and without advection) showed no variation in peak water level, and no damping was exhibited.
The fundamental seiche period for a rectangular basin is given by:
|
(C2) |
where l is the length of the basin. For a wave traveling parallel to the x-axis of the slosh-test grid, T = 0.02337 day, corresponding to a frequency of 42.8 cycles/day (cpd). Spectrums of water levels shown in Figures C3 and C4 for Tests 1 and 2 are displayed in Figures C5 and C6, respectively. Peak amplitudes occur at 42.8 cpd. Small peaks (slightly visible in Figures C5 and C6) are contained in the spectra at 85.6 cpd, indicating the presence of a weak harmonic of the primary seiche frequency. These tests indicate that the discretized equations applied for calculation of water level in Tests 1 and 2 do not introduce erroneous frequencies into the solution.
Figure C SEQ AppAFig \* MERGEFORMAT 3. Time series of water level for Test 1
Figure C SEQ AppAFig \* MERGEFORMAT 4. Time series of water level for Test 2
Figure C SEQ AppAFig \* MERGEFORMAT 5. Spectrum of water level for Test 1
Figure C SEQ AppAFig \* MERGEFORMAT 6. Spectrum of water level for Test 2
Tests 3 and 4 were designed to determine how damping and mass conservation properties respond to a grid with sloped bottom and constant cell sizes. Grid configurations and numerical parameters for Tests 3 and 4 are given in Table C3 and the grid is shown in Figure C7. Changes in volume of 5 ´ 10-6 and 1 ´ 10-5 percent were calculated for Tests 3 and 4, respectively.
Table C Slosh Test Grid Configurations and Numerical Parameters for Tests 3 and 4 | ||
Numerical Parameter | Value | |
D'''s''', m | 500 | |
D'''t''', s | 5 | |
h, 'm' | Min = 9 m, Max = 11 m, Dh /Dx = 0.0002 | |
Length, m | 10,000 | |
Width, m | 2,500 |
Figure C SEQ AppAFig \* MERGEFORMAT 7. Computational grid for Tests 3 and 4
Damping properties for Tests 3 and 4 are shown by time series of water level in Figures C8 and C9, respectively. Starting maximum water level for both tests was 0.01 m at 0.5Dx from the end of the grid. Test 3 does not exhibit damping over the duration of the simulation. Water levels from Test 4 are slightly damped by properties of the numerical algorithm for the advective terms. In both tests, water levels are weakly nonlinear, as described for Tests 1 and 2.
Spectrums of water levels for Tests 3 and 4 are displayed in Figures C10 and C11, respectively. Peak amplitudes occur at 42.4 cpd. Small peaks (slightly visible in Figures C5 and C6) are contained in the spectra at 84.8 cpd, indicating the presence of a weak harmonic of the primary seiche frequency. These spectra demonstrate that the discretized equations applied for calculation of water level in Tests 3 and 4 do not introduce erroneous frequencies into the solution.
Figure C SEQ AppAFig \* MERGEFORMAT 8. Time series of water level for Test 3
Figure C SEQ AppAFig \* MERGEFORMAT 9. Time series of water level for Test 4
Figure C SEQ AppAFig \* MERGEFORMAT 10. Spectrum of water level for Test 3
Figure C SEQ AppAFig \* MERGEFORMAT 11. Spectrum of water level for Test 4
Appendix D
Wind Set Up Tests
This appendix describes tests designed to evaluate model performance under wind forcing. Two tests were performed to examine model response. The first test compares the results of numerical calculations to those of an analytical solution for a rectangular channel. The second test evaluates model response to wind forcing in a basin of arbitrary shape.
Wind set up in a rectangular channel for steady-state conditions is maintained by a balance of the pressure gradient, wind forcing, and bottom friction. For 1D flow, the momentum equation (Equation 2) can be reduced to (Bretschneider 1966)[2]:
|
(D1)
|
where
is the wind stress, and
is the bottom stress. In the absence of bottom friction, Equation D1 of Bretschneider can be rewritten as:
5+ 9
|
(D2)
|
where
. The coefficient k was derived from wind and bottom-stress coefficients determined in a study of Lake Okeechobee, FL (Saville 1953). Application of Equation D2 to the present study required modification of k to be equal to the wind-stress coefficient applied in CMS-M2D, which varies with wind speed, so that analytical and numerical solutions could be compared. Integration of Equation D2 gives:
|
(D3)
|
where C is a constant of integration that is determined by the position of the channel origin if no bottom is exposed, and by the position of the water-land interface if drying has occurred. The parabolic form of Equation D3 indicates a locally wind-forced water surface in a rectangular channel does not have a linear distribution along the channel axis. Instead, the water surface will be curved, with the curvature being dependent on the channel length and depth, and the wind speed.
For comparison of the analytical and numerical solutions of Equation D2, a 1D grid was developed that had a total length of 9,500 m and was 95 cells long. Each square cell had side dimensions of 100 m. The still-water depth was set to 2 m and the time step was 10 s. The wind speed set to 10 m/s and held constant. A 10-m/s wind speed corresponds to
. Figure D1 compares water level computed analytically from Equation D2 and numerically by CMS-M2D, under steady-state conditions. The numerical solution well reproduces the steady-state, wind-induced set up and set down of the analytical solution for simulations applied to shallow water.
Figure D1. | Comparison of analytical and numerical solutions of steady-state, wind-induced set up, no exposed channel bottom
|
To test the response of the model to wind forcing in a basin of arbitrary geometry, a modified version of the hourglass test developed for a curvilinear grid (Spaulding 1984) was examined. An hourglass-shaped grid was developed with constant cell-side lengths set to 500 m and an ambient depth of 5 m. The time step was set to 30 s. Two tests were performed, one for wind blowing from the west (from the negative x-axis) and the other for wind blowing from the north (from the positive y-axis). For each test, the wind was ramped up over a 24-hr interval and held constant at 10 m/s for the remainder of the 100-hr simulation.
Water-surface elevation for the west wind test is shown in Figure D2. The set up and set down are symmetric between the upper and lower portions of the grid. The water-surface elevation contours for the north wind test are shown in Figure D3. The water-surface elevation for the north wind test is equally distributed between the left and right portions of the grid. Symmetry in both tests demonstrates that there is no inconsistency in calculations at wall boundaries. The results of this test demonstrate that the model responds to wind forcing with qualitatively correct patterns of water-surface deviation. Wind set up calculated by CMS-M2D is qualitatively similar to results displayed graphically by Spaulding (1984) and by Cialone (1989).
Figure D2. Water-surface elevation contours for west wind
Figure D3. Water-surface elevation contours for north wind
Appendix E
Advection Test
Advective terms were examined by developing a grid with an idealized jettied inlet and displaying the current field during ebb flow. Inclusion of the advective terms in the momentum equations allows the model to properly simulate eddies and an ebb jet as the ebb flow exits the inlet. Eddies and jets will not develop without advective terms represented in the calculations.
The idealized inlet grid was specified to have a flat bottom with a still-water depth of 5 m. Inlet dimensions were 300 m wide and 800 m long. Cell-side lengths were 25 m on all sides in the inlet and nearshore areas. Cell size increased to 50 m offshore. The time step was set to 1 s. Water level was forced with a sine wave having amplitude of 0.5 m and period of 12.42 hr.
Two simulations were conducted, one without advective terms included in the computations and one with them. Figure E1 shows the current field from the test without the advective terms during the ebb tidal cycle. No eddies are present and flow exiting the inlet spreads laterally, but does not form a jet. However, with the advective terms included in the computations (Figure E2), a well-developed ebb jet is formed and eddies are present on each side of the jet. The ebb jet and eddies are symmetrical about the axis of the inlet, demonstrating that there is no directional bias in the advection algorithm.
Figure E1. Ebb current field for advection test, no advective terms
Figure E2. Ebb current field for advection test, with advective terms
Appendix F
Wetting and Drying Tests
Tests were performed on a sloped-bottom grid to determine the stability properties of the wetting and drying algorithms. One- and two-dimensional tests were conducted to demonstrate the behavior of the flooding and drying algorithms under different sources of hydrodynamic forcing and basin geometries. A list and short description of the tests are given in Table F1. Comprehensive descriptions and model results are provided for specific tests.
Table F1 List and Description of Flooding and Drying Tests | ||
Test Number | Title | Description |
1 | Channel with sloping bottom | Test wetting and drying of a 1D canal forced by time-varying wind. |
2 | Depression in a square grid | Test stability for wetting dry cells surrounding a ponded area. 2D grid. |
Test 1: Channel with Sloping Bottom
A 1D test was conducted to examine the flooding and drying algorithms for a channel with simple geometry. The test consisted of a channel with a sloping bottom that was forced by a time-varying wind. Channel length was 150,000 m with 101 cells along the main channel axis (parallel to the x-axis) and 1 grid-cell wide. Cells were square with side lengths of 1,500 m. The time step was 100 s. Channel bottom slope was
with maximum and minimum still-water depths of 3 m and 0.5 m, respectively. Forcing consisted of a wind that blew along the channel from the shallow end toward the deep end. Wind was ramped up over 10 hr to a maximum speed of 10 m/s and held constant from hour 10 to hour 100. Wind was turned off at 101 hr and remained off for the remainder of the 200-hr simulation.
Time series of water elevation for the sloping-bottom test are shown in Figure F1 in 10-hr intervals. Twenty-six cells on the shallow (left) end of the grid are shown out of 101 cells in the grid. Dry cells are shown as white. When wind stress is applied to the channel (0 to 100 hr), water is forced out of the shallow area so that successive drying of cells occurs from left to right. A maximum of 21 cells dried during the simulation. After the wind was turned off (hour 101), water flowed from right to left, successively inundating the dry cells.
Figure F SEQ AppFFig \* MERGEFORMAT \* MERGEFORMAT 1. Demonstration of flooding and drying of sloped channel
Test 2: Depression in Square Grid
A source of instability in some flooding and drying algorithms is rewetting of dry cells that surround an area that has formed a pond. A robust flooding algorithm will remain stable with any pattern of wet and dry cells. This test was designed to determine the behavior of the wetting and drying algorithms where ponding takes place.
A square grid consisting of 400 50 ´ 50-m cells was developed in which the central area contained a square ring of shallow cells surrounding deeper cells (Figure F2). Outside the shallow ring, the depth was set to 2 m everywhere. Cells in the shallow ring were specified with depths of 0.3 m, with adjacent interior cells having depths of 0.5 m, and adjacent outer cells having depths of 1.0 m. Four cells in the center of the depression were 1-m deep.
Figure F2. Grid with central depression surrounded by shallow area
The grid was forced with water-surface elevation at six cells located on the left side of the grid, denoted by the red line in Figure F2. A sinusoidal forcing was applied having amplitude of 1 m and period of 12 hr. The depth at which cells are deemed dry was specified as 0.05 m. The model was run for duration of 50 hr, with a ramp of 1 day and time step of 5 s. Plots of water level and circulation are shown in Figures F3 through F10 at hourly intervals starting at hour 42. During hours 43 to 47, the water level was lowered to expose the ring of shallow cells. The dry cells are denoted by white with no current vectors. A pond was formed within the interior of the ring of dry cells.
Water level rises from hours 45 through 49. At hour 46 (Figure F7), cells on the outer portion of the shallow ring (with ambient depths of 1 m) have rewetted. By hour 48 (Figure F9), all previously dry cells have rewetted. During the process of flooding and drying, no instabilities were encountered. Thus, the wetting and drying algorithms in CMS-M2D function properly for the situation of flooding around a ponded area.
Figure F3. Water level and circulation pattern at hour 42
Figure F4. Water level and circulation pattern at hour 43
Figure F5. Water level and circulation pattern at hour 44
Figure F6. Water level and circulation pattern at hour 45
Figure F7. Water level and circulation pattern at hour 46
Figure F8. Water level and circulation pattern at hour 47
Figure F9. Water level and circulation pattern at hour 48
Figure F10. Water level and circulation pattern at hour 49
Appendix G
Wave-Adjusted Boundary Condition Test
This appendix demonstrates the wave-adjusted boundary condition for an idealized beach. In the example, wave forcing is applied to calculate a longshore current. The wave-adjusted boundary condition modifies the water level and velocity on boundaries forced by water-surface elevation to account for the presence of wave forcing. Thus, the prescribed boundary forcing is applied after being adjusted for the presence of waves. This adjustment promotes realistic circulation patterns and wave set up and set down at the boundaries, making values at the lateral boundaries compatible with values calculated within the calculation domain near the boundaries. A simulation without the wave-adjusted boundary condition is also shown, which illustrates the necessity of the modification of water level and velocity at the boundary.
An idealized beach computational grid was developed (Figure G1) having grid spacing of 50 m in each horizontal direction, and depths ranging from 1.13 m at the shoreline to 18.4 m at the seaward boundary. CMS-M2D and STWAVE domains were identical for the idealized beach simulations.
Forcing for the models was specified as:
a. | STWAVE: Wave height of 2.0 m, period of 10.0 s, and wave angle of 25.0 deg. Spectral spreading parameters were
|
b. | CMS-M2D: Water-surface elevation set to 0.0 m along lateral and seaward boundaries.
|
CMS-M2D was run for 12.5 hr with a 2-s time step and ramp duration of 0.02 day. Results from the simulation without the wave-adjusted boundary condition are described first, followed by those with the wave-adjusted boundary condition.
Figure G1. Computational grid for idealized beach
Water-surface elevation and velocity as computed by the standard water-surface elevation forcing boundary condition (boundaries not adjusted for presence of wave forcing) are shown in Figures G2 and G3, respectively. Set up is calculated at the shoreline, but decreases to zero at the lateral boundaries. A longshore current is develops in the nearshore area, but velocities at the lateral boundaries are weaker than inside the calculation domain and display unrealistic patterns. The unrealistic velocity patterns and water-surface elevations at the lateral boundaries owe to the absence of the wave stresses there.
Water-surface elevation and velocity as computed by the wave-adjusted boundary condition are shown in Figures G4 and G5, respectively. Set up is calculated at the shoreline and extends nearly uniformly across the entire grid. An alongshore current is well developed in the nearshore area, both near the lateral boundaries and away from them. Because the boundary conditions were adjusted for the presence of wave forcing, the circulation patterns are realistic. Water can flow onto and off of the grid along the lateral boundaries without incompatibility between the boundary condition and interior hydrodynamic processes.
Figure G2. | Water-surface elevation calculated with standard water surface elevation boundary condition
|
Figure G3. | Velocity calculated with standard water-surface elevation boundary condition
|
Figure G4. | Water-surface elevation calculated with wave-adjusted boundary condition
|
Figure G5. Velocity calculated with wave-adjusted boundary condition
Appendix H
Surf Zone Verification Tests
To verify CMS-M2D’s response to wave forcing, calculations were compared to laboratory data collected by Visser (1982).[3] Because the laboratory experiments were conducted with monochromatic waves and the bathymetry was uniform alongshore, CMS-M2D was modified to include a 1D monochromatic wave-transformation model. This model applies linear wave theory to compute wave transformation and computes wave height of depth-limited breaking and broken waves by the wave stability relationship of Dally et al. (1985). It also includes representation of the wave roller as described by Dally and Brown (1995) and Dally and Osiecki (1994). Water levels calculated by CMS-M2D are passed to the wave transformation model to allow waves to respond to the set up and set down.
Larson and Kraus (2002) provide succinct descriptions of the wave stability relationship and roller model together with comparisons of calculated and measured longshore current, including comparisons to the Visser (1982) data, with and without the roller model. Their results demonstrate that inclusion of the roller momentum for calculation of the longshore current and wave set up increases accuracy. Without the roller momentum, peak longshore current was calculated to be seaward of the measured peak, whereas including the roller shifted the peak toward the beach and aligned it with the measured peak. This adjustment of the cross-shore distribution of the longshore current modified the set up with an increased set up with the roller as compared to without it.
Calculation of longshore current and set up are conducted for two of the Visser (1982) data sets, Cases 4 and 7. Wave properties for each case are listed in Table H1, and cross-shore profiles and roughness information are listed in Table H2. The major difference between these two cases is the bottom roughness. Case 4 was run on a smooth concrete bottom, whereas bottom roughness elements were introduced for Case 7, with the wave conditions held to be the same as much as possible. By choosing these two tests, performance of CMS-M2D in the surf zone over a great difference in bottom roughness can be demonstrated.
Table H1 Wave Properties at Offshore End of Basin and Breaker-to-Depth Ratio for Selected Visser (1982) Laboratory Experiments | ||||
Case | Wave Height, m | Wave Period, s | Wave Angle, deg | Breaker to Depth Ratio |
4 | 0.076 | 1.02 | 19.9 | 0.83 |
7 | 0.078 | 1.02 | 18.4 | 0.74
|
Table H2 Cross-Shore Profiles and Roughness Information for Selected Visser (1982) Laboratory Experiments | ||||
Case | Distance Across Shore, m | Depth, m | Beach Slope1 | Bottom Roughness |
4 | -0.472 7.328 8.253 |
-0.040 0.350 0.350 |
0.050 | smooth |
7 | -0.472 7.328 9.503 |
-0.040 0.350 0.350 |
0.050 | rough |
1 Beach slope calculated for sloping portion of basin.
|
Computational grids were developed for each bathymetric profile, and depths were uniform alongshore. For each grid, the cell spacing parallel to the alongshore direction (
) was 0.10 m. Cell spacing parallel to the cross-shore direction (
) was 0.01 m from the offshore boundary to approximately 1 m from the shoreline. Resolution was increased nearer to the shore to capture details of the set up calculation. As an example, resolution change for the Case 7 computational grid is shown in Figure H1.
Simulations were conducted with identical parameter settings and time step, with the exception of the Manning roughness coefficient and the breaker-to-depth ratio for which settings were case dependent. The depth at which cells were specified as dry was 0.001 m and the time step was 0.025 s. Breaker-to-depth ratios for each case are provided in Table H1. The Manning roughness coefficient, serving as an adjustment parameter, was set to 0.0125 for Case 4 and to 0.022 for Case 7, after numerical experimentation. Wave calculations proceeded shoreward until a minimum wave height of 0.001 m was reached.
Figure H1. | Detail of computational grid for Case 7 showing cross-shore variation in cell spacing
|
Comparison of measurements and calculations for wave height, longshore current speed, and set up for Case 4 are shown in Figures H2, H3, and H4, respectively. Calculated wave height compared well with measurements, except for the maximum wave height, which was underestimated by the wave-transformation model. The calculated breakpoint took place directly shoreward and nearly at the measured location.
Calculated longshore current speed well reproduced the distribution determined by the measurements (Figure H3). Peak current speeds compare well, with the calculated peak located directly seaward of the measured peak. Seaward of the breakpoint, calculated weakening of the longshore current closely follows the measured reduction in current.
Set up calculations also compare favorably to the measurements (Figure H4). The slope of the set up is slightly reduced in the calculations as compared to measurements, but the overall result indicates that the model is reproducing the wave-induced set up.
Comparison of calculated and measured wave height and longshore current speed for Case 7 (artificial roughness) are shown in Figures H5 and H6, respectively, and calculated set up is shown in Figure H7 (set up measurements are not available for Case 7). Calculated peak wave height underpredicts the measured peak wave height. However, because measured wave heights before breaking are not available, there is no indication of whether strong nonlinear shoaling took place. Calculated longshore current speed closely matches measurements over the peak and offshore reduction in speed.
Figure H2. Wave height comparison for Case 4
Figure H3. Longshore current comparison for Case 4
Figure H4. Wave-induced set up comparison for Case 4
Figure H5. Wave height calculated by roller model for Case 7
Figure H6. Longshore current calculated by M2D for Case 7
Figure H7. Wave-induced set up calculated by M2D for Case 7
In conclusion, calculations with CMS-M2D well reproduced the qualitative and quantitative features of accurate measurements of the longshore current in the laboratory. This test validates the overall equation solution scheme of the model and contributions of several terms, including the radiation stresses, bottom friction, and mixing.
Appendix I
Hard-Bottom Tests
Functioning of the hard-bottom algorithm within CMS-M2D is demonstrated through two idealized tests that are designed to show calculation properties of the capability. Specifically, demonstrated properties are:
a. | Hard-bottom cells do not erode below the depth of the non-erodible substrate.
|
b. | Symmetry in depth change calculations for symmetric flow and bathymetry.
|
c. | Hard-bottom cells can accrete and erode sand above the non-erodible substrate.
|
d. | Sand volume is conserved.
|
Both tests applied the Watanabe total load sediment transport formulation, and sand size was specified as 0.2 mm.
Test 1: Hard-Bottom Calculations Under Symmetric Conditions
Test 1 was designed to demonstrate that the hard-bottom algorithm does not allow erosion of hard substrate, but does allow accumulation of sand over the non-erodible surface. This example also demonstrates symmetry of the methodology and capability to represent complex configurations of non-erodible cells. The bathymetric configuration (Figure I1) was specified to be a plane horizontal bottom with a pyramid-shaped shoal in the center. The configuration was subject to sediment transport generated by a constant and uniform current in four different simulations where the current originated from each of the four sides of the square. All simulations gave the same symmetrical result relative to the main current direction, so that results from only one direction are provided.
Figure I1. Model bathymetry with cells shown
Figure I2 shows the simulation where the current is directed from the left side. Because of the presence of the pyramidal sediment mound, the current bends around the structure as indicated in the figure. Under the influence of the current, the pyramid erodes.
Figure I3 shows a magnification of the center of the grid. Cells with triangles are specified as having a non-erodible bottom at the initial cell depth (i.e., no sediment cover on the hard-bottom layer). Two sets of runs are compared. In the first set, the hard bottom is not in place, i.e., allowing the cells to erode. In the second set, the hard-bottom cells are in place, limiting erosion to a specified hard-bottom depth. Erosion below the hard-bottom level is allowed in these cells during the first set of runs, but not during the second set.
The result of this comparison is shown in Figure I4, displaying the change in depth (blue denotes erosion, and yellow/red denote accretion) after 7 hr of simulation with no hard-bottom cells present (triangles are retained in the image for reference to previous figure). Cells specified as hard bottom in the following set of simulations erode in the present simulation, indicating that the hydrodynamic conditions cause the bottom to erode below the level of the hard-bottom depth specified in the following simulation. The depth change pattern is symmetric. Figure I5 shows the same situation, but with the hard bottom in place. The calculations show that there is no erosion below the hard-bottom cells. Also, symmetry is preserved.
Figure I2. Velocity field over the schematized bottom
Figure I3. | Magnification of bathymetry at the center of grid. Cells with triangles are specified as having hard bottom
|
Although not shown here, other simulations with the current originating from the other three sides of the square gave the same symmetrical results relative to the main current direction. In all simulations, sand volume was conserved.
Figure I4. | Change in depth after 7 hr of simulation without hard-bottom cells (triangles are left in image for reference to previous figure). Blue denotes erosion; yellow/red denote accretion
|
Figure I5. | Change in depth after 7 hr of simulation with hard-bottom cells. Blue denotes erosion and yellow/red denote accretion
|
Test 2: Erosion and Accretion of Sand Over Hard Substrate
A computational grid was developed representing a canal with a shallow region located in its center (Figure I6). Grid depths vary along the x-axis and are uniform along the y-axis. In Figure I6, cells having hard substrates are denoted with triangles. During the simulation, these cells will erode to the hard-bottom surface, then undergo deposition and erosion. The hard- bottom depth was specified to be 0.7 m, and initial grid depths at these cells were 0.5 m, meaning that 0.2 m of sand covered the hard substrate at the beginning of the simulation.
Figure I6. | Computational grid and bathymetry for demonstration of erosion and deposition of sand over hard substrate. Triangles denote hard-bottom cells
|
A 400-hr simulation was conducted in which the current flowed to the right (positive) for 200 hr and then to the left (negative) for 200 hr, establishing conditions for both erosion and deposition at the hard-bottom cells. Time series of current velocity and depth at one hard-bottom cell are shown in Figure I7. During the first 90 hr of the simulation, scour of the sand layer takes place. Once the hard-bottom substrate is reached, scour ceases and the depth remains constant at 0.7 m until current velocities are weak enough for shoaling to start (approximately hour 210). Sand accumulates until the current is strong enough to begin erosion (approximately hour 280), and the erosion continues until the hard substrate is reached again. Thus, the hard-bottom algorithm functions correctly for situations in which sand overlays the hard substrate and erosion takes place to the level of the hard bottom, and for situations in which conditions become depositional at hard-bottom cells.
Figure I7. | Time series of depth of bottom and velocity for demonstration of scour and accumulation of sand over hard substrate
|
- ↑ References cited in this appendix are included in the References section at the end of the main text of this report.
- ↑ All references cited in this appendix are included in the References section at the end of the main text.
- ↑ All references cited in this appendix are included in the References section at the end of the main text.