Sediment Model

Sediment Model

Overview

The ROMS model represents sediment using separate cohesive and non-cohesive categories, each with an unlimited number of user-defined size classes. Each class has fixed attributes of grain diameter, density, settling velocity, critical shear stress for erosion, and erodibility constant. These properties are used to help determine the bulk properties of each bed layer. Most of the information in this section is also presented in Warner, Sherwood, Signell, Harris, and Arango, 2008 "Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model", Computers & Geosciences.

Cohesive sediment dynamics (i.e. consolidation, flocculation) have been implemented in development branches and are presently being migrated to the trunk. For more information, contact Chris Sherwood.

Sediment bed

The sediment bed is represented by three-dimensional arrays with a fixed number of layers beneath each horizontal model cell. Each cell of each layer in the bed is initialized with a thickness, sediment-class distribution, porosity, and age. The mass of each sediment class in each cell can be determined from these values and the grain density. The bed framework also includes two-dimensional arrays that describe the evolving properties of the seabed, including bulk properties of the surface layer (active layer thickness, mean grain diameter, mean density, mean settling velocity, mean critical stress for erosion) and descriptions of the subgrid scale morphology (ripple height and wavelength). These properties are used to estimate bed roughness in the BBL formulations and feed into the bottom stress calculations. The bottom stresses are then used by the sediment routines to determine resuspension and transport, providing a feedback from the sediment dynamics to the hydrodynamics.

The bed layers are modified at each time step to account for erosion and deposition and track stratigraphy. At the beginning of each time step, an active layer thickness za is calculated based on the relation of Harris and Wiberg (1997). The thickness of the top bed layer has a minimum thickness equivalent to za. If the top layer is thicker than za, no action is required. If the top layer is less than za thick, then the top layer thickness is increased by entraining sediment mass from deeper layers until the top layer thickness equals za. If sediment from deeper than the second layer is mixed into the top layer, the bottom layer is split to enforce a constant number of layers and conservation of sediment mass. Each sediment class can be transported by suspended-load and/or bedload (described below). Suspended-load mass is exchanged vertically between the water column and the top bed layer. Mass of each sediment class available for transport is limited to the mass available in the active layer. Bedload mass is exchanged horizontally between the top layers of the bed. Mass of each sediment class available for transport is limited to the mass available in the top layer. Suspended-sediment that is deposited, or bedload that is transported into a computational cell, is added to the top bed layer. If continuous deposition results in a top layer thicker than a user-defined threshold, a new layer is provided to begin accumulation of depositing mass. The bottom two layers are then combined to conserve the number of layers. After erosion and deposition have been calculated, the active-layer thickness is recalculated and bed layers readjusted to accommodate it. This step mixes away any very thin layer (less than the active layer thickness) of newly deposited material. Finally the surficial sediment characteristics, such as D50, ripple geometry, etc., are updated and made available to the bottom stress calculations.

Suspended-sediment transport

Sediment suspended in the water column is transported, like other conservative tracers (e.g., temperature and salinity) by solving the advection–diffusion equation (5) with a source/sink term for vertical settling and erosion. The vertical advection algorithm includes a piece-wise parabolic method (Colella and Woodward, 1984) and a weighted essentially non-oscillatory (WENO) scheme (Liu et al., 1994). This method allows the integration bounds of depositional flux to use multiple grid boxes in the vertical direction, so it is not constrained by the CFL criterion. Zero-flux boundary conditions are imposed at the surface and bottom in the vertical diffusion equation. The source or sink term in the advection equation represents the net of upward flux of eroded material and downward settling.

Bedload transport

This version of ROMS implements two methods for computing bedload transport: 1) the Meyer-Peter Müeller (1948) formulation for unidirectional flow and 2) the formulae of Soulsby and Damgaard (2005) that accounts for combined effects of currents and waves. The formulae depend on the characteristics of individual sediment classes, including median size $D_{50}$, grain density $\rho _{s}$, specific density in water $s=\rho _{s}/\rho$, and critical shear stress $\tau _{c}$. Non-dimensional transport rates $\Phi$ are calculated for each sediment class and converted to dimensional bedload transport rates $q_{bl}$ using

$$q_{bl}=\Phi {\sqrt {\left({s-1}\right)gD_{50}^{3}}}\rho _{s}$$

These are horizontal vector quantities with directions that correspond to the combined bed-stress vectors.

Details of each of these bedload-transport methods are discussed here:

• Meyer-Peter Müeller
• Soulsby & Damgaard

Morphology

The bed model accounts for changes in sea floor elevation resulting from convergence or divergence in sediment fluxes. These morphological changes can have significant influence on flow and transport when they are larger than a few percent of the water depth. The morphological changes are accounted for by equating the bottom boundary condition of the vertical velocity to the rate of change of elevation of the sea floor. This method is completely mass conserving and retains tracer constancy preservation.

A morphological scale factor is also provided to allow an increased rate of morphological change, which can be useful for simulating evolution over long time periods. Strategies for morphological updating are described by Roelvink (2006). In our implementation, bedload fluxes, erosion, and deposition rates are multiplied by a scale factor. A scale factor with a value of one has no effect, and values greater than one accelerate the bed response. For bedload transport, the scale factor is multiplied against the bedload transport rates. For suspended load transport, the scale factor multiplies the exchange of sediment (erosive or depositional flux) at the bed-water interface. The magnitude of sediment concentrations in the water column are not modified – just the exchange rate to and from the bed. For both bedload and suspended load, sediment is limited in availability as described previously, based on the true amount of sediment mass (not multiplied by the scale factor). This morphological scale factor method works well for systems with unlimited sediment in the bed. However, it can generate extra sediment in systems with limited supplies of bed sediment. This occurs when the amount of sediment to be eroded is limited by the amount available and application of the morphological scale factor cannot remove the scaled amount of sediment from the bed. Subsequent deposition does place a scaled amount of sediment on the bed thus creating new mass in the bed. Other approach (Lesser et al 2004) is to limit the amount of sediment fluxed to the water column in these situations. This gives unrealistically low sediment concentrations, but conserves bed sediment.

Sediment Density

Effects of suspended sediment on the density field are included with terms for the weight of each sediment class in the equation of state for seawater density as

$$\rho =\rho _{water}+\sum \nolimits _{m=1}^{Nsed}{{\frac {C_{m}}{\rho _{s,m}}}\left({\rho _{s,m}-\rho _{water}}\right)}$$

This enables the model to simulate processes where sediment density influences hydrodynamics, such as density stratification and gravitationally driven flows.

Bottom Stress

Reynolds stresses, production and dissipation of turbulent kinetic energy, and gradients in velocity and suspended-sediment concentrations vary over short vertical distances, especially near the bed, and can be difficult to resolve with the vertical grid spacing used in regional-scale applications. ROMS provides algorithms to parameterize some of these subgrid-scale processes in the water column and in the bottom boundary layer (BBL). Treatment of the BBL is important for the circulation model solution because it determines the stress exerted on the flow by the bottom, which enters the Reynolds-averaged Navier-Stokes equations as a boundary conditions for momentum in the $x$ and $y$ directions:

$$K_{M}{\frac {\partial u}{\partial s}}=\tau _{bx}^{};\,\,\,\,\,\,K_{M}{\frac {\partial v}{\partial s}}=\tau _{by}^{}$$

Determination of the BBL is even more important for the sediment-transport formulations because bottom stress determines the transport rate for bedload and the resuspension rate for suspended sediment.

ROMS implements either of two methods for representing BBL processes: (1) simple drag-coefficient expressions or (2) more complex formulations that represent the interactions of wave and currents over a moveable bed. The drag-coefficient methods implement formulae for linear bottom friction, quadratic bottom friction, or a logarithmic profile. The other, more complex methods, implement some of the many wave-current BBL models (e.g., Jonsson and Carlsen, 1976; Smith, 1977; Grant and Madsen, 1979; Madsen, 1994; Styles and Glenn, 2000) and couple them with calculations of bottom roughness. ROMS offers three methods that implement slightly different combinations of algorithms for the wave-current interactions and moveable bed roughness. The first method (SG_BBL) is based on the wave-current algorithm and the ripple geometry and moveable bed roughness of Styles and Glenn (2000, 2002). The second method (MB_BBL) uses efficient wave-current BBL computations developed by Soulsby (1995) in combination with sediment and bedform roughness estimates of Grant and Madsen (1982), Nielsen (1986) and Li and Amos (2001). These algorithms and an example of their use on the Southern California continental shelf are described by Blaas et al. (2005). The third method (SSW_BBL) implements either the wave-current BBL model of Madsen (1994) or that of Styles and Glenn (2000) along with moveable bed routines proposed by Wiberg and Harris (1994; Harris and Wiberg, 2001). The differences in approach among these routines are small, but they can produce significantly different results.

The linear and quadratic drag-coefficient methods depend only on velocity components $u$ and $v$ in the bottom grid cell and constant, spatially-uniform coefficients $\gamma _{1}$ and $\gamma _{2}$ specified as input:

$$\tau _{bx}^{}=\left({\gamma _{1}+\gamma _{2}{\sqrt {u^{2}+v^{2}}}}\right)u$$

$$\tau _{by}^{}=\left({\gamma _{1}+\gamma _{2}{\sqrt {u^{2}+v^{2}}}}\right)v$$

where $\gamma _{1}$ is the linear drag coefficient and $\gamma _{2}$ is the quadratic drag coefficient. The user can choose between linear or quadratic drag by setting one of these coefficients to zero. The bottom stresses computed from these formulae depend on the elevation of $u$ and $v$ (computed at the vertical mid-elevation of the bottom computational cell). Therefore, in this s-coordinate model, the same drag coefficient will be imposed throughout the domain even though the vertical location of the velocity is different.

The logarithmic formulation assumes that flow in the BBL has the classic vertical logarithmic profile defined by a shear velocity $u_{*}$ and bottom roughness length $z_{0}$ as

$$|u|={\frac {u_{*}}{\kappa }}\ln \left({\frac {z}{z_{0}}}\right)$$

where speed $|u|={\sqrt {u^{2}+v^{2}}}$ , friction velocity $u_{*}={\sqrt {\left|{\tau _{bx}^{}}\right|+\left|{\tau _{by}^{}}\right|}}$, $z$ is the elevation above the bottom (vertical mid-elevation point of the bottom cell), $\kappa =0.41$ is von Kármán’s constant, and $z_{0}$ is a constant (but possibly spatially varying) bottom roughness length (m). Kinematic stresses are calculated as:

$$\tau _{bx}^{}={\frac {\kappa ^{2}u{\sqrt {u^{2}+v^{2}}}}{\ln ^{2}(z/z_{0})}}$$

$$\tau _{by}^{}={\frac {\kappa ^{2}v{\sqrt {u^{2}+v^{2}}}}{\ln ^{2}(z/z_{0})}}$$

The advantage of this approach is that the velocity and the vertical elevation of that velocity are used in the equation. Because the vertical elevation of the velocity in the bottom computational cell will vary spatially and temporally, the inclusion of the elevation provides a more consistent formulation.

More complex routines are required to simulate BBL processes in the presence of waves and mobile sediment. The short (order 10-s) oscillatory shear of wave-induced motions in a thin (a few cm) wave-boundary layer produces turbulence and generates large instantaneous shear stresses. The turbulence enhances momentum transfer, effectively increasing the coupling between the flow and the bottom and increasing the frictional drag exerted on the mean flow, averaged over many wave periods. The large instantaneous shear stresses often dominate sediment resuspension and enhance bedload transport. Sediment transport can remold the bed into ripples and other bedforms, which present roughness elements to the flow. Bedload transport can also induce drag on the flow, because momentum is transferred to particles as they are removed from the bed and accelerated by the flow. Resuspended sediments can cause sediment-induced stratification and, at high concentrations, change the effective viscosity of the fluid.

The BBL parameterization implemented in ROMS requires inputs of velocities $u$ and $v$ at reference elevation $z_{r}$, representative wave-orbital velocity amplitude $u_{b}$, wave period $T$, and wave propagation direction $\theta$ (degrees, clockwise from north). The wave parameters may be the output of a wave model such as SWAN or simpler calculations based on specified surface wave parameters and should represent the full spectrum of motion near the bed (cf. Madsen, 1994; Wiberg and Sherwood, 2008). Additionally the BBL models require bottom sediment characteristics (median grain diameter $D_{50}$, mean sediment density $\rho _{s}$, and representative settling velocity $w_{s}$); these are based on the composition of the uppermost active layer of the bed sediment during the previous time step. Bed stresses associated with mean current above the wave-boundary layer $\tau _{b}$, wave motions $\tau _{w}$, and maximum vector sum of the two $\tau _{w}c$ from the previous time step are used as initial estimates for the next time level.

Descriptions of alternative bottom-boundary layer formulations are described here.

• SSW_BBL
• MB_BBL
• SG_BBL