See below for extensive information on the following:

- Elliptic Grid Generator
- Hydrodynamic Model
- Helical Flow
- Sediment Transport
- Bed Load
- Suspended Load
- Bed Forms and Hydraulic Resistance
- Bank Erosion
- Large Scale Morphology
- Other Modules

A description of the parameters to be specified in the MIKE 11 ST and MIKE 21 C input file menus is also available. You can download the entire help system for MIKE 21 C in the compiled HTML (*.chm) format.

**ELLIPTIC GRID GENERATOR**

The standard version of the MIKE 21 is based on a rectilinear computational grid. For modeling of the open sea, and most estuaries and coastal applications, such a grid provides sufficient accuracy. In river application, however, an accurate resolution of the boundaries is required which necessitates for the use of curvilinear or non-structured grids. Curvilinear grids have the advantage over the non-structured grids that the computational schemes are much faster.

The benefits of using a curvilinear computational grid compared to a rectilinear grid are visualized in Figure 1, where a river reach is described in both a curvilinear and a rectilinear grid. In this example, the discretization in the curvilinear grid is based on 210 computational points, whereas the rectilinear model uses 228 active (water) points. The curvilinear model provides a much better resolution of the flow near the boundaries and thereby a higher modeling accuracy. Bigger time step can be used in the curvilinear model because grid lines follow the streamlines. In the rectilinear model, a total of 900 points is defined and stored, whereas for the curvilinear model, only 270 points are defined. Thus, the required storage capacity is relatively larger in a rectilinear model.

MIKE 21 C is based on a so-called orthogonal curvilinear grid. This is created with a grid generator, which solves a set of elliptic partial differential equations. The advantage of using an orthogonal grid is that the finite difference equations describing the two-dimensional flow become substantially simpler than if a general (non-orthogonal) curvilinear grid was applied. This implies that the numerical scheme becomes more accurate with an orthogonal grid and that the computational speed of the engine improves.

The orthogonal curvilinear grid used in MIKE 21 C is obtained from the following elliptic partial differential equations:

where

- x,y Cartesian co-ordinates
- s,n curvilinear co-ordinates (anti-clockwise system)
- g "weight" function

The weight function is a measure of the ratio between the grid cell length in the s- and the n-direction, respectively. It is defined through the following relationships:

The boundary condition for this system is the non-linear orthogonality condition:

where the former expresses the condition of orthogonality and the latter expresses the location of grid points (x,y) on a specific curve describing the boundary. A strongly implicit method is employed for solving the partial differential equations with a special Newton-Raphson procedure for the boundary conditions, see Stone (1968). The input for the grid generator is the grid weights, g, in every grid point and the boundary conditions in terms of a set of (xb,yb) co-ordinates for each of four boundary lines.

Output from the grid generator comprises the (x,y) co-ordinates of the grid line intersection points. The weight function, g, can be linked to the water depth or similar hydraulic parameters in such a way that the grid line density depends on the local hydraulic conditions. Figure 2 shows an example of a curvilinear grid obtained with the elliptic grid generator.

Generation of an orthogonal curvilinear grid is an iterative process in which boundaries are smoothened and weight functions are adjusted until the computational grid is judged to be sound, i.e. without too large gradients in grid cell spacing and curvature of grid lines. Initial hydrodynamic simulations may still reveal the need for further refinements and adjustments of the grid. The developed graphical facility is an important tool for quick plotting of the computational grid and in assessing the quality of the generated grids.

Figure 2 Example of curvilinear grid. The weight function causes the grid to be concentrated in the deeper parts (depth adaptive grid) (click to enlarge).

**HYDRODYNAMIC MODEL**

The hydrodynamic model simulates the water level variation and flows in rivers and estuaries. Model simulations are based on a curvilinear computational grid covering the area of interest. The MIKE 21 hydrodynamic model has been under continuous development at the Danish Hydraulic Institute since 1970 whereas the curvilinear version of MIKE 21 using the same solution methods was developed in 1990.

The hydrodynamic model solves the full dynamic and vertically integrated equations of continuity and conservation of momentum (the Saint Venant equations) in two directions. The following effects can be included in the equations when used for river applications:

- flow acceleration
- convective and cross-momentum
- pressure gradients (water surface slopes)
- bed shear stress
- momentum dispersion (through e.g. the Smagorinsky formulation)
- Coriolis forces
- wind forces
- flow curvature and helical flow

The curvature of the grid lines gives rise to additional terms in the partial differential equations for the flow. The equations solved in MIKE 21 C are:

where

- s,n co-ordinates in the curvilinear co-ordinate system
- p,q mass fluxes in the s- and n-direction, respectively
- H water level
- h water depth
- g gravitational acceleration
- C Chezy roughness coefficient
- Rs, Rn radius of curvature of s- and n-line, respectively.
- RHS The right hand side describing a.o. Reynold stresses, Coriolis force, wind friction, atmospheric pressure (see MIKE 21 Hydrodynamic Module)

Additional terms due to the curvilinear co-ordinate system also emerge when describing the Reynold stresses.

Figure 3 Space staggered grid used in the hydrodynamic model of MIKE 21 C (click to enlarge)

The setting up of a hydrodynamic river model using MIKE 21 C involves the following three steps:

- The extent of the modeling area is selected and the computational grid is designed for the river. If a curvilinear grid is applied, the grid generation described in the previous section is used.
- The bathymetry is entered into the computational grid, i.e. bed levels are specified in each computational grid point. Various utility programs in MIKE 21 PP performing conversions, interpolation, extrapolation, smoothening etc. are available.
- The specification of the boundary conditions for instance a time series of discharge at the inflow boundary, and a time series of water level at the downstream model boundary.

The calibration process for the MIKE 21 C hydrodynamic model involves tuning of a number of calibration factors. All the calibration factors have a physical meaning and should not be arbitrarily given values outside their realistic ranges to obtain agreement with observed data.

**HELICAL FLOW**

Helical flow is a principal secondary flow phenomenon in rivers. Whilst it does not have any strong influence on the general flow pattern in rivers with large width-depth ratios, it has a significant influence on the sediment transport direction and hence the morphological changes in the river channel (see Olesen, 1987). The helical flow is, therefore, calculated in connection with sediment transport simulations when larger scale morphology is modeled. It is an important ingredient in the development of bend scour, confluence scour, and in formation of point bars as well as alternating bars.

The helical flows occur in curved flows, especially in river bends. It arises from the imbalance between the pressure gradient and the centripetal acceleration working on a water particle moving along a curved path. Near the riverbed the helical flow is directed towards the center of flow curvature. The magnitude of the helical flow (i.e. the transverse flow velocity component) rarely exceeds 5-10 % of the main flow velocity in natural rivers. A typical helical flow pattern is shown in Figure 6.

The helical flows occur in curved flows, especially in river bends. It arises from the imbalance between the pressure gradient and the centripetal acceleration working on a water particle moving along a curved path. Near the riverbed the helical flow is directed towards the center of flow curvature. The magnitude of the helical flow (i.e. the transverse flow velocity component) rarely exceeds 5-10 % of the main flow velocity in natural rivers.

Assuming a logarithmic main flow velocity distribution over the vertical and a parabolic eddy viscosity distribution, the magnitude of the secondary flow can be shown to be proportional to the main flow velocity, the depth of flow and the curvature of the main flow stream lines. In MIKE 21 C, the streamline curvature is calculated explicitly from the depth integrated flow field. The gradual adaptation of secondary (helical) flow to changing curvature is accounted for by solving a first order differential equation along the stream lines with the strength of the helical flow as the dependent variable (see also de Vriend, 1981). The strength of the helical flow is used to determine the direction of both bed and suspended load.

With streamlines of radius Rs of curvature and depth H, the helical flow intensity can be expressed as (U is the main flow speed, C is Chezys number):

where the helical flow strength was first derived by Rozowsky, 1957:

**SEDIMENT TRANSPORT**

In connection with detailed (two-dimensional) mathematical modeling of sediment transport and morphology in rivers with large suspended load transports, it is necessary to distinguish between bed and suspended load in order to:

- simulate the dynamic development of bed form dimensions
- account for the effect of helical flow as well as the bed slope on the sediment transport direction

The relatively simple total load sediment transport formulas, such as Engelund & Hansen, 1967 and Ackers & White, 1973 can therefore not be used for river applications with consideration of helical flow and bed slope without a separate specification of how the sediment is distributed between bed load and suspended load. It is possible, though, to run simulations by MIKE 21 C with a total transport formula only by disregarding the effect of helical flow and bed slope.

The sediment transport models developed by Engelund & Fredsoe, 1976, and van Rijn, 1984, which distinguish between bed and suspended load form the basis for the sediment transport description in MIKE 21 C. The specification of the sediment transport formulas is, however, very flexible. Specially developed sediment transport formulas (determined from for instance field measurements) can be specified separately for bed load and suspended load, respectively. With this flexible sediment transport formulation, it is also possible to select formulas like Engelund-Hansen, Smart-Jaeggi, and Meyer-Peter.

The MIKE 21 C is capable of simulating several fractions of sediment, each with a characteristic grain size. Each grain fraction is simulated separately with due consideration to the interaction between the various grain components at the riverbed (shielding, armoring). With this model, it is possible to simulate the sorting in space and time of graded sediment, which is essential in many river applications, especially in gravel river models.

If the suspended sediment transport is negligible compared to the bed load transport, the suspended sediment model can be switched off, so only a bed load model (or a total load model) is employed.

**BED LOAD**

In MIKE 21 C, the bed load transport is calculated explicitly from one of the selected formulas, e.g. Engelund-Fredsøe, van Rijn, or Meyer-Peter. These formulas relate the transport rate to the bed shear stress and the grain diameter.

On a horizontal bed, the transport direction will coincide with the direction of the bed shear stress. The direction of the bed shear stress may, however, deviate from the direction of the depth-averaged flow due to the helical flow as described in a previous section. If d is the angle of deviation due to the helical flow and c is the helical flow strength, the relation tan d = c exists.

The bed slope is I, the angle of deviation is a , q is non-dimensional bed shear stress and G and a are calibration parameters. Typical values of G and a are 0.66 and 0.5, respectively. Most of these relations have, however, only been verified against data from laboratory tests and are, therefore, not necessarily applicable to natural rivers with large suspended load, see also Talmon/Struiksma/van Mierlo 1995. The applicable relation is therefore often determined via calibration of the model. Corrections in the calculated bed load transport due to sloping river bed in the main flow direction is also done in MIKE 21 C.

Figure 7 Shear stress on bed load sediment particles is composed of three "drag" components: Force due to main flow Fbs, force due to secondary flow Fbn, force due to gravity on a sloping riverbed (click to enlarge).

The output from the bed load model is the bed load transport rate and direction at every computational grid point.

**SUSPENDED LOAD**

Standard methods for calculation of suspended load are not applicable in the case of detailed (ie. high resolution) modeling of rivers. It is necessary to include the time-space lag in the sediment transport response to changes in local hydraulic conditions. Consider, for instance, an increase of flow velocity in a river (e.g. a constriction). As the flow velocity increases, the entrainment at the river bed will increase correspondingly, but it will take some time (and hence some distance) for the sediment entrained at the bottom to disperse all through the water column. This means that the actual suspended load is not only a function of local hydraulic conditions, as it is normally assumed in most mathematical sediment transport models, but it is also a function of what takes place upstream and earlier in time.

A relevant time scale for the time lag (T) is the settling time for a sediment grain in the water column:

T = h*/ws

where h* is the effective fall height (depends on the shape of the vertical concentration profile and, thus, on the fall velocity and eddy dispersion) and ws is the fall velocity.

Correspondingly, a length scale for the space lag is

L = T× U

where U is the flow velocity. Assuming h*=8 m, ws=0.02 m/s, U =2 m/s, the time and length scales are 400 s and 800 m, respectively. Consequently, the space lag is important for river applications .

The space lag effect is modeled in MIKE 21 C by a depth-averaged convection-dispersion model which represents the transport and the vertical distribution of suspended solids and flow. This model is an extension of the model first developed for one dimension by Galappatti (1983) and later in two dimensions by Wang (1989). The model by Wang, however, was primarily developed for estuarine applications and did not include the effect of helical flow.

The secondary flow profile is computed in MIKE 21 C and used together with the primary flow profile and the concentration profile when the suspended load is integrated over the depth. As the concentration is highest near the river bed, the suspended load transport will be deflected towards the center of flow curvature. In contrast to bed load, the transverse river bed slope does not influence the direction of suspended load transport.

The depth averaged convection-dispersion model requires an expression for the equilibrium concentration. The models by eg. Engelund & Fredsøe (1982) or van Rijn (1984) can be used for that purpose. The empirical formulas implemented in MIKE 21 C can also be used assuming that the equilibrium concentration equals the suspended load transport divided by the water flux.

The output from the suspended load model is the concentration of suspended sediment as well as the suspended load transport and direction at every computational grid point.

**BED FORMS AND HYDRAULICS RESISTANCE**

It is far more complex to determine the hydraulic resistance in alluvial rivers than in channels with a fixed bed. This is because a large part of the hydraulic resistance in alluvial rivers is caused by bed form drag. The bed forms have a configuration determined by the sediment transport and flow. The hydraulic resistance will therefore exhibit both temporal and spatial variations.

Generally, the hydraulic resistance is divided into that caused by drag on the bed forms (form friction) and due to shear forces on the bed (skin friction). The skin friction is relatively accurately determined from a logarithmic boundary layer equation based on the median grain size of the bed. The form friction, however, can only be determined analytically if the size of the bed forms is known.

Several hydraulic resistance predictors for alluvial rivers have been proposed, including those of Engelund-Hansen and Ackers & White. Both of these models are semi-empirical linking the hydraulic resistance to the local instantaneous flow conditions. However, there can be a significant lag between the form friction (i.e. the bed form size) and the hydraulic conditions. For instance, in many tropical rivers, with distinct dry and wet periods, the sand dunes found at the river bed during the dry season have been formed during a receding flood at the end of the flood season. In such rivers, the hydraulic resistance will often be relatively larger during the dry season. It is therefore important to account for this time lag.

In MIKE 21 C, the dynamic development of the bed form size (height and length) is calculated with the model developed by Fredsøe (1979). In a subsequent step, the form friction is calculated using a Carnot type formula for expansion loss. The skin friction is determined from a logarithmic boundary layer equation. In quasi-steady flow, this model suggests that as the flow velocity increases, the dune size and hence the hydraulic resistance increases. As the flow velocity increases further, the bed form height and water depth first increase at more or less the same rate and hence the hydraulic resistance only changes slowly until the sand dunes are washed away relatively abruptly and the hydraulic resistance decreases rapidly.

For some applications, the differences in bed resistance due to rapidly changing bed levels may be more important than the differences in bed resistance due to increase of decrease in the flow rate. Thus, use of a simpler alluvial bed resistance model was found to be beneficial for such cases. This model reads (C is the Chezy number, h is the local depth, a and b are calibration constants):

C = a hb

For b=1/6, the model simply equals the Manning formula, whereas for b=0, the model is the Chezy formula. In the morphological model, the depth h will change both as a result of bed level changes and due to water level changes.

**BANK EROSION**

An important aspect of river morphological processes is bank erosion. For natural rivers without protection of the riverbanks, appropriate description of bank erosion may be vital to get a full picture of the overall morphology of the area. In MIKE 21 C, the bank erosion can be simulated in parallel with the sediment transport and hydrodynamic simulations. In each time step, the eroded bank material is included in the solution of the sediment continuity equation.

The following bank erosion model is incorporated:

where

- Eb the bank erosion rate in m/s
- z local bed level
- S near bank sediment transport
- h local water depth
- a,b,g calibration coefficients specified in the model

The extra sediment which is discharged into the river from this source and which is included in the sediment continuity equation is (hb height of the bank above the water level):

The accumulated bank erosion causes a retreat of the bank line position and thereby also in the extent of the modeling area. The MIKE 21 C incorporates these plan form changes by re-generating the curvilinear grid during the simulation when the bank line changes exceed a certain pre-defined threshold. In this manner, the morphological model becomes a plan form model.

For an example of a simulation with updating of the curvilinear grid, see details of the tutorial of a meandering stream.

**LARGE SCALE MORPHOLOGY**

Large scale in this context means changes in the general bed level within the model boundaries as well as plan form changes due to bank erosion.

The bed and suspended load sediment transport models described earlier in this section predict sediment transport rate and direction. The change of bed level is, therefore, easily determined by integrating the net inflow or outflow within a control volume causing either deposition or scour.

This integration is explicit, due to the non-linear character of the sediment transport relations applied. This implies that the large-scale morphological model needs to be subjected to a rigorous stability criterion for the time step. However, a time step substantially larger than the time step in the hydrodynamic model can generally be applied. In addition, the MIKE 21 C has the capability to simulate hydrodynamics assuming quasi-steadiness. With this feature, it is possible to simulate very long time series (years) with even very detailed 2-D models. Finally, the time step for sediment transport simulation can vary automatically by specifying a maximum acceptable sediment transport courant number rather than the time step it self. This option is useful when simulating time series with large variations from low to high flow conditions.

Three kinds of boundary conditions can be specified for the morphological model at the upstream boundary: Total sediment transport, bed level changes, or concentration of suspended sediment (in which case, bed load is calculated explicitly from local hydraulic condition at the boundary). The boundary conditions may vary both in space and time along the boundaries. The effect of bank erosion is included in the sediment continuity equation, ie. the eroded bank material is included as a lateral boundary condition.

The output from the large-scale morphological model is sediment transport (bed load as well as suspended load), the bed level and bed level changes at every grid point and at every time step. Also accumulated bank erosion and new grid co-ordinates are output from the model if the bank erosion and grid update module is activated.

**OTHER MODULES**

The MIKE 21 C presently contains the following modules

- Hydrodynamics, HD
- Advection-Dispersion, AD
- Sediment transport in rivers, ST
- Mud transport, MT

The HD, AD and ST modules are described here. The various water quality modules known from MIKE 21 can be installed with little effort as the core of these modules is the AD model, which is already available in curvilinear co-ordinates. No wave modules are available in MIKE 21 C, nor the specific ST modules required for simulating coastal processes (e.g. interaction of waves and currents).

The MT module in MIKE 21 C is similar to the standard MIKE 21 MT module, which describes the erosion, transport and deposition of cohesive sediments (mud, silt and clay) under the action of waves and currents (wave modeling to be carried out separately with MIKE 21 NSW). The model also takes into account the consolidation of the riverbed. The advection-dispersion scheme used in the AD-module is used to describe the transport and dispersion of suspended sediments. The model can be used to determine the siltation of cohesive materials in harbors, lagoons, rivers and estuaries and to determine the fate of dredged spoils.