Stabilized Formulation for Modeling the Erosion/Deposition Flux of Sediment in Circulation/CFD Models

Abstract

In field-scale modeling, when the resuspension of sediment is modeled using a hydrodynamic model, a standard and common approach is to add a resuspension flux as the bottom boundary condition in the transport model. In this study, we show that the way of simply imposing an empirical bottom erosion formula as the flux is actually unrealistic. Its inability to stabilize the sediment concentration can cause excessive suspension fluxes in some extreme cases. Moreover, we present a modified erosion/deposition formula to model the resuspension of sediment. The formulation is based on volume conservation in the presence of erosion/deposition near the bottom. By taking into account the prescribed dry density of the bed material, the proposed formulation is able to produce realistic near-bed concentrations while ensuring model stability. The formulation is then tested in a one-dimensional vertical model and field modeling cases using a three-dimensional coastal circulation model. We show that the modified formulation is particularly important in modeling mud resuspension subject to the large bottom stress, which can be a result of waves or a strong river discharge.

1. Introduction

Sediment transport is an important factor influencing estuarine and coastal environments. Dynamics of sediment budgets tightly interact with hydrological process, morphology, and ecological functions in a wide range of spatial and temporal scales. A proper sediment transport model is valuable in guiding the management of water resources in coastal and estuarine areas while facing alterations due to human activities and climate change. The modeling of suspended sediment transport typically involves a hydrodynamics model and sediment transport model. While hydrodynamics modeling has remarkably progressed recently due to the advancing computation power and techniques, modeling sediment transport remains challenging. From the modeling perspective, sediment processes can simply be divided into two regimes [1]. The first regime is found in the regions with dense concentrations of particles, in which particle-particle interactions are intense, resulting in the dominance of friction and collision. The second regime is found in areas of dilute concentration such that hydrodynamic forcing dominates particle motion. The transport of sediment in the first regime is referred to as bed load while the second regime is suspended load. To deal with sophisticated particle-particle interactions, modeling bed load transport relies on empirical parameterization (e.g., [2]). When particle size is small, such as fine sand, silt, and clay, modeling suspended sediment is analogous to the transport of scalar concentration field, C, in conjunction with settling velocity, $w_s$ [3].

Unlike bed load, which is usually connected to morphodynamics, contexts of suspended load are usually more related to hydrodynamic problems in stratified fluid, such as flow instabilities and turbulence suppression. Therefore, typical field modeling for suspended sediment transport usually disregards the associated morphodynamic change (e.g., [4,5]). The typical governing equation describing the evolution of suspended sediment concentration (SSC) is an advection-diffusion equation, which is written as [3].

(1)$\frac{\partial C}{\partial t}+\nabla \times \left[\left(U-w_s{\widehat{e}}_3\right)C\right]=K_T{\nabla }^{2}C,$

where $U=\left(u,v,w\right)$ is the velocity field obtained from hydrodynamics calculation, and $K_T$ is the eddy viscosity. Equation (1) is a standard approach widely employed in conjunction with various field-scale hydrodynamics models to describe suspended sediment transport. When coupled with the hydrodynamics model, additional sink or source terms in Equation (1) are required to account for deposition or erosion, respectively. Thus, disregarding any lateral transport and vertical turbulent diffusion, an integral form of mass conservation within the bottommost cell is written as

(2)${\int }_{c.v.}\frac{{\displaystyle \partial C}}{{\displaystyle \partial t}}dV={\oint }_{c.s.}{F}_{b}ds,$

where $c.v.$ indicates the control volume, $c.s.$ indicates the surface of the control volume, and $F_b>0$ represents erosion and <0 represents deposition. For a grid cell with vertical size $\Delta z$, a sink/source term to be applied to the bottommost cell can be obtained by reducing Equation (2) to

(3)$\frac{{\displaystyle \partial \overline{C}}}{{\displaystyle \partial t}}=\frac{{\displaystyle F_b}}{{\displaystyle \Delta z}},$

where the overbar indicates the volume average. With a validated hydrodynamics model, the predictability of Equation (1) relies on the appropriate empirical formulas to parameterize sediment inflow flux from the bottom. A formula used to describe the inflow flux from the bottom is called the pick-up function in the context of sand transport [6] or erosion formula in the context of mud resuspension [7]. It is usually formulated based on laboratory flume studies (e.g., [8,9,10,11]) or in situ flume experiment (e.g., [10,12,13,14,15]) as a function of the difference between the bottom shear stress and the critical shear stress of erosion. The erosion rate is usually described as a power-law [15,16,17] or exponential [8,9,10] relationship with the stress difference for the hard bed or the soft bed erosion, respectively. Up to now, the formulation that has been most widely used is the power law with the number of power = 1, which gives the linear relation between the erosion and stress difference, based on the laboratory results of Ariathurai and Arulanandan [18]. It has been used for suspended sediment transport modeling in conjunction with a number of renowned hydrodynamics models, such as the Regional Ocean Modeling System (ROMS) [4,19,20,21,22,23,24,25,26,27], Unstructured Tidal, Residual, Intertidal Mudflat (UnTRIM) model [28], and Finite Volume Community Ocean Model (FVCOM) [29,30]. It disregards the erosion of the top fluffy mud layer (i.e., soft bed erosion). Moreover, several experimental studies show that the power number for the hard bed erosion ranges from 2 to 6. That is, for the more realistic modeling for the bottom resuspension of sediment, an erosion formula that is more sensitive to the stress difference (i.e., either the exponential law or power law with higher power number) needs to be adopted, e.g., the mud transport module in MIKE21 [31]. However, no matter what types of the erosion formula is employed, implementation must be undertaken with caution. A lack of constraints preventing excessive sediment erosion can result in over-resuspension in the model, particular when the erosion formula is sensitive to the stress difference. This is of critical importance in situations involving wave action or strong river discharge, which can lead to strong bed shear stress and erosion. In such cases, excessive erosion that is overpredicted by the model can result in erroneous predictions of sediment transport.

In fact, the simple addition of empirical sink and source terms to Equation (1) is not representative of the actual cases. It should be noted that the erosion and deposition are also associated with the changes in the depth of the water column. Although this change can be small compared to the total depth, as demonstrated in the present study, it plays a key role in maintaining the magnitude of near-bed SSC below that of the dry density of the bed material. This enables predictions of greater precision and stabilizes the simulations. This study reexamines the mechanisms underlying the erosion and deposition while taking into account the conservation of volume on a moving-grid framework. Our derivation shows that the moving-grid calculations provide a connection to the dry density of the bed. In the case of erosion, this prevents excessive erosion and maintains near-bed SSC at values below those of the prescribed dry density of the bed. Based on the moving-grid formulation, we then present modified sink/source terms for regular fixed-grid calculation. Through numerical testing, we demonstrate that the proposed modified sink/source terms are able to produce realistic near-bed SSC. Without a lengthy discussion related to the physical details of sediment transport, the present study emphasizes the importance of considering changes in the depth of the water column during erosion and deposition. It also presents a modified formulation that guarantees model stability while giving realistic values of SSC.

The rest of this article is organized as follows. In Section 2, mathematical expressions of erosion and deposition are derived, and modified sink and source terms are presented. In Section 3, we perform numerical tests using both moving and fixed grids. Moreover, an example of three-dimensional modeling for San Francisco Bay, USA, is presented. The summary and conclusion are then given in Section 4.

2. Modified Erosion/Deposition Equations Based on Moving-Bottom Formulation

Existing modeling frameworks with fixed grids often lead to discrepancies between the model and actual conditions, in which the depth of water slightly increases with erosion and decreases with deposition. To demonstrate this, we start from the general formulation of the Reynolds transport theorem to describe conservation of sediment mass within a control volume (denoted with $c.v.$):

(4)$\frac{{\displaystyle d}}{{\displaystyle dt}}{\int }_{c.v.}CdV+{\oint }_{c.s.}C{u}_{rel}\times dA=0,$

where the second integral on the left-hand side (LHS) indicates that the integral is taken with respect to the surface of the control volume, and $u_{rel}$ is the velocity relative to the control surface. Erosion and deposition on the bed lead to a corresponding change in control volume. In the following, we deal with cases of erosion and deposition separately.

2.1. Deposition

Neglecting flow velocities that advect sediment, Equation (4) describing the settling of sediment onto the bed can be rewritten as

(5)$\frac{{\displaystyle d}}{{\displaystyle dt}}{\int }_{c.v.}CdV+{\oint }_{c.s.}C\left(-w_s-u_g\right)\times dA=0,$

where $u_g$ is the moving velocity of the control surface due to deposition. As depicted in Figure 1b, disregarding any vertical flux at the top boundary of the control volume, temporal discretization of Equation (5) using the explicit Euler scheme for a control volume of height $\Delta z$ leads to

(6)$C^{n+1}(\Delta z-h_d)-C^n\Delta z+C_b^n\left(w_s+u_{g,b}\right)\Delta t=0,$

where C is the volume-averaged concentration, $C_b$ is the concentration at the bottom boundary, the superscript n represents the computational time step, $h_d$ is the (positive) height of deposition from time = t to $t+\Delta t$, and $u_{g,b}$ is the (positive) moving speed of the bottom of the control volume. Treating the bottommost grid cell in the computation domain as the control volume of Equation (6), the concentration within the bottommost cell, denoted as $C_1$, can be obtained by rearranging Equation (6) as

(7)$C_1^{n+1}=C_1^n-C_b^n\frac{{\displaystyle w_s\Delta t}}{\Delta z}\frac{{\displaystyle \Delta z}}{\Delta z-h_d}-\left(C_b^n-C_1^n\right)\frac{h_d}{\Delta z-h_d}$

The third term on the LHS of Equation (6) is the total mass of deposition. Given that the bed material has a dry density, $C_{dry}$, the height of deposition, $h_d$, is thereby expressed as $C_{dry}{h}_d=C_b\left(w_s+u_g\right)\Delta t$. Using the identity, $u_{g,b}\Delta t=h_d$, $h_d$ can be obtained with

(8)$h_d=\frac{{\displaystyle C_b^n}}{{\displaystyle C_{dry}-C_b^n}}{w}_s\Delta t$

Equation (7) describes the “realistic” volume-averaged concentration resulting from only deposition in the bottommost cell, within which conservation of volume is regarded. If $h_d$ is sufficiently small to be disregarded, Equation (7) is reduced to the standard discretized equation on fixed grids to obtain the concentration in the bottommost cell, denoted using the overhat, as

(9)${\widehat{C}}_1^{n+1}=C_1^n-C_b^n\frac{{\displaystyle w_s\Delta t}}{\Delta z}$

Substitution of Equation (9) into Equation (7) leads to

(10)$C_1^{n+1}={\widehat{C}}_1^{n+1}-\frac{{\displaystyle h_d}}{{\displaystyle \Delta z-h_d}}\left[C_b^n\left(1+\frac{{\displaystyle w_s\Delta t}}{\Delta z}\right)-C_1^n\right]$

In numerical calculations, $C_b$ is usually extrapolated from interior values, and, in most cases, $C_b>C_1$. Therefore, Equation (10) indicates that, in the case of pure deposition, as long as the deposition height does not exceed height of the cell, the realistic concentration must be less than the calculated value obtained using the standard deposition formula (Equation (9), as depicted in Figure 1a).

2.2. Erosion

The case of erosion differs slightly from that of deposition. One can see this by starting from the general formulation of mass conservation for pure erosion, which is expressed as

(11)$\frac{{\displaystyle d}}{{\displaystyle dt}}{\int }_{c.v.}CdV+{\oint }_{c.s.}C\left(u^{*}-u_g\right)\times dA=0,$

where $u^{*}$ parameterizes a fluctuation velocity that results in the suspension flux of bed sediment into the water column. Due to the fact that complex processes involved in bottom erosion cannot be resolved in any field-scale hydrodynamics model, the second term on the LHS of Equation (11) relies on empirical formulas based on laboratory or field measurements. Practically, the measured mass flux, as denoted by $E_b$, represents the magnitude of the second integral of Equation (11) (i.e., the net mass flux). Using the same time discretization as in Section 2.1 and disregarding any flux on the top of the control volume, as depicted in Figure 1d, one obtains $C^{n+1}(\Delta z+h_e)-C^n\Delta z-E_b\Delta t=0$, where $h_e$ is the (positive) erosion depth and can be obtained with

(12)$h_e=\frac{{\displaystyle E_b\Delta t}}{{\displaystyle C_{dry}}}$

Again, the volume-averaged concentration in the bottommost cell (denoted as $C_1$) can be obtained with

(13)$C_1^{n+1}=C_1^n\frac{{\displaystyle \Delta z}}{{\displaystyle \Delta z+h_e}}+\frac{{\displaystyle E_b\Delta t}}{\Delta z+h_e}$

It can be seen that, if $h_e$ is negligible, Equation (13) is reduced to the standard erosion equation describing changes in concentration due to erosion in the bottommost grid cell, i.e.,

(14)${\widehat{C}}_1^{n+1}=C_1^n+\frac{{\displaystyle E_b\Delta t}}{\Delta z}$

Substitution of Equation (14) into Equation (13) leads to

(15)$C_1^{n+1}={\widehat{C}}_1^{n+1}-\frac{{\displaystyle h_e}}{{\displaystyle \Delta z+h_e}}\left(\frac{{\displaystyle E_b\Delta t}}{{\displaystyle \Delta z}}+C_1^n\right),$

which shows that, with the occurrence of erosion, the resulting near-bottom concentration is always less than the value obtained using the standard erosion formulation (Equation (14)). Although $h_e$ can be sufficiently small so that the difference between Equations (14) and (15) is negligible, the most important feature of Equation (15) is the fact that considering the depth of erosion provides a linkage between the near-bottom concentration ($C_1$) and the dry density of the bed ($C_{dry}$). One can see this by rearranging Equation (15) to become

(16)$C_1^{n+1}=C_{dry}\left[\frac{{\displaystyle C_1^n\Delta z+E_b\Delta t}}{{\displaystyle C_{dry}\Delta z+E_b\Delta t}}\right]$

Therefore, as $C_{dry}$ is always higher than $C_1$, Equation (16) shows that concentration in the bottommost cell would never exceed the concentration of the bed material. However, this cannot be guaranteed in the standard fixed-grid implementation (Equation (14)), as depicted in Figure 1c, in which the connection to bed properties is missing. As a result, the boundary condition of bottom erosion is ill-posed in the sense that it can continuously provide sediment through the bottom boundary regardless of the near-bed concentration. This situation can be far worse when power-law or exponential-law empirical formulas are adopted for the parameterization of erosion due to the nonlinear growth of the eroded quantity with the increasing excessive bottom shear stress (e.g., Equation (21)).

In the conditions of strong turbulent mixing, such as strong tidal flows, the resulting stability issue may not be severe because eroded sediments can be soon suspended into the water column and vertically well mixed, but it would not be so when wave forcing is involved. The bottom shear stress resulting from the wave-induced oscillating boundary layer is usually much stronger than those induced by tidal currents only. Moreover, in shallow water estuaries, the wave effect is relatively strong when water depth is low, such as low tides, during which the turbulence intensity is weak. As a result, once eroded, sediments can accumulate in the bottommost cell and exceed the dry density of the bed. When this occurs, the resulting strong buoyant flux of sediment can result in numerical instability and inaccuracy.

2.3. Modified Erosion/Deposition Formula

Based on the previous derivation, we propose the “modified flux” (denoted with the superscript $\ast$) to reformulate the bottom sink/source term (Equation (3)) on fixed grids as

(17)$\frac{{\displaystyle \partial \overline{C}}}{{\displaystyle \partial t}}=\frac{{\displaystyle F_b^{*}}}{{\displaystyle \Delta z}},$

in which $F_b^{*}$ is derived in accordance with changes in the control volume (Equations (7) and (13)), to obtain

(18)$F_b^{*}=\left\{\begin{array}{cc}-C_b^n{w}_s\left[\frac{{\displaystyle \Delta z}}{{\displaystyle \Delta z-h_d}}\left(\frac{{\displaystyle C_{dry}-C_1^n}}{{\displaystyle C_{dry}-C_b^n}}\right)\right]\hfill & ,\tau <{\tau }_{CD}\left(\mathrm{deposition}\right),\\ E_b\left[\frac{{\displaystyle \Delta z}}{{\displaystyle \Delta z+h_e}}\left(1-\frac{{\displaystyle C_1^n}}{{\displaystyle C_{dry}}}\right)\right]\hfill & ,\tau >{\tau }_{CE}\left(\mathrm{erosion}\right),\end{array}\right.$

where $\tau$ is the bed shear stress based on hydrodynamic calculation, ${\tau }_{CD}$ is the critical shear stress for deposition, and ${\tau }_{CE}$ ($\ge {\tau }_{CD}$) is the critical shear stress for erosion. In Equation (18), $h_d$ and $h_e$ are calculated using Equations (8) and (12), respectively. It should be noted that Equation (17) does not conserve mass at the water-bottom interface but rather maintains a stable, realistic concentration in near-bed regions. Therefore, if the calculation of deposited or eroded mass is necessary in the cases requiring the estimation of the erosion/deposition depth, the actual eroded mass flux is $E_b$ and deposited mass flux is $w_s{C}_b$.

3. Numerical Tests

3.1. 1DV Model

In this section, we conduct simple numerical tests showing the differences among standard erosion formulation, the moving-grid case, and the modified erosion formulation. Because the present study focuses on vertical suspension and deposition, we solve one-dimensional vertical (1DV) equations disregarding any lateral transport. In conjunction with a turbulence closure model, the resulting formulation provides a standard approach to obtain the vertical distribution of SSC in horizontally homogeneous turbulent flow. When taking into account changes in water depth due to erosion and deposition, the 1DV model is solved on the moving frame, with a moving-grid velocity, $w_g$. Thus, the equations for horizontal momentum (u) and sediment concentration (C) are respectively written as

(19)$\frac{{\displaystyle \partial u}}{{\displaystyle \partial t}}-\frac{{\displaystyle \partial }}{{\displaystyle \partial z}}\left(w_{g}u\right)=-\frac{{\displaystyle 1}}{{\displaystyle \rho }}\frac{{\displaystyle \partial p}}{{\displaystyle \partial x}}+\frac{{\displaystyle \partial }}{{\displaystyle \partial z}}\left[\left(\nu +{\nu }_T\right)\frac{{\displaystyle \partial u}}{{\displaystyle \partial z}}\right]$

and

(20)$\frac{{\displaystyle \partial C}}{{\displaystyle \partial t}}-\frac{{\displaystyle \partial }}{{\displaystyle \partial z}}\left(w_{s}C\right)-\frac{{\displaystyle \partial }}{{\displaystyle \partial z}}\left(w_{g}C\right)=\frac{{\displaystyle \partial }}{{\displaystyle \partial z}}\left[\frac{{\displaystyle {\nu }_T}}{{\displaystyle {\sigma }_T}}\frac{{\displaystyle \partial C}}{{\displaystyle \partial z}}\right],$

where $\nu$ is the kinematic viscosity of water, ${\nu }_T$ is the eddy viscosity, p is the pressure, x is the horizontal coordinate, and ${\sigma }_T=0.7$ is the turbulent Schmidt number. The $k-ϵ$ model is employed for the turbulence closure that also takes the buoyancy effect due to suspended sediments into account [32]. In Equation (19) through Equation (20), $w_g$ is non-zero only in the moving-grid calculation. In that case, the domain depth varies at each time step according to Equations (8) and (12), and new equal-spacing grid cells are obtained accordingly using the same total grid number. Using finite-volume discretization in space, the resulting numerical scheme is an arbitrary Lagrangian Eulerian (ALE) scheme capable of guaranteeing mass conservation [33], representing a simplified version of three-dimensional flow modeling coupled with the moving bed form in Chou and Fringer [34]. As previously mentioned, source and sink terms are added to the bottommost cell of sediment concentration calculation. Here, we use the formulation proposed by Parchure and Mehta [8] for erosion, which is expressed as

(21)$E_b=E_{b,0}exp\left[\alpha {\left(\tau -{\tau }_{CE}\right)}^{\beta }\right],$

where $E_{b,0}$ is the erodibility, $\tau$ is the bottom shear stress, ${\tau }_{CE}$ is the critical shear stress of erosion, and $\alpha$ and $\beta$ are two site-specific constant coefficients. For deposition, we set ${\tau }_{CD}={\tau }_{CE}$ in Equation (18), i.e., sediments begin to deposit where there is no erosion. According to the field case of cohesive sediment modeling by Lumborg and Pejrup [35] and the values suggested in MIKE 21 [31], this study uses $E_0=0.01$ g m⁻² s⁻¹, $\alpha =6.5$, $\beta =1.0$, and ${\tau }_{CE}=0.3$ N m⁻², which represents erosion over the soft mud layer. Although the choice of these values is arbitrary, it does not affect the conclusions of the present study. Just beneath the water column is a partially-consolidated mud layer with $C_{dry}=530$ g L⁻¹. In Equation (21), $\tau$ is the result of combined wave current forces. Assuming that the wave orbital velocity is aligned with the current velocity, $\tau$ can be written as [36]

(22)$\tau =\sqrt{{\tau }_c^2+{\tau }_w^2+2{\tau }_c{\tau }_{w}cos{\theta }_{cw}},$

where ${\tau }_c$ is the shear stress induced by currents, ${\tau }_w$ is induced by waves, and ${\theta }_{cw}$ is the angle between currents and waves and is set to 1 here in the 1DV model. In coastal circulation models, ${\tau }_c$ can be obtained using the log law with a prescribed value for bottom roughness ($z_0$), and ${\tau }_w$ is calculated based on wave properties that can be obtained through the coupled wave model (e.g., [19]). Thus, ${\tau }_w$ is written as

(23)${\tau }_w=\frac{{\displaystyle 1}}{{\displaystyle 2}}\rho f_w{\left|u_{orb}\right|}^2,$

where $f_w$ is the wave friction coefficient and $u_{orb}$ is the wave orbital velocity.

In many cases, ${\tau }_w$ can be much higher than ${\tau }_c$ such that waves become the major contributor to sediment erosion. However, further suspension into the water column depends on turbulence induced by tidal currents. Therefore, it can be expected that the most unstable case would occur during strong wave events when tidal currents are weak (e.g., during low-water slack tide). Here, for the sake of simplicity, we assigned arbitrary values directly to ${\tau }_w$ rather than considering wave properties.

Numerical simulations are performed in a one-dimensional vertical domain with depth of 3 m, discretized using 20 equally-spaced grid cells. A settling velocity of $w_s=2\times {10}^{-4}$ m s⁻¹ is used corresponding to the particle diameter $d_0\approx 46$ µm based on the Stokes law. A steady current is forced by a pressure gradient resulting in a depth-averaged steady current velocity = 0.2 m s⁻¹. Using the log law with $z_0=d_0/12$, the simulated ${\tau }_c$ is equal to $0.0014$ N m⁻², which is less than ${\tau }_{CE}$ and would not induce any suspension. Three cases associated with different wave conditions are simulated using $u_{orb}=0.55$, 0.70, and 0.82 m s⁻¹. Assuming that the wave period is 10 s and the wavelength is 30 m, thtese $u_{orb}$s correspond to wave heights of 1.17, 1.49, and 1.75 m, respectively. Using the formulation for $f_w$ given in Madsen et al. [37], these three cases correspond to ${\tau }_w=0.85$, 1.33, and 1.74 N m⁻². For each wave condition, simulations are conducted according to various modeling approaches. The first is the moving-grid simulation using the original deposition/erosion formulation (Equation (2)). Because it strictly conserves volume, it gives the most realistic and accurate SSC results. Both the second and third modeling approaches are performed on fixed grids. In the second case, we apply the modified deposition/erosion formulation (Equations (17) and (18)). In the third case, the standard approach that directly employs the original deposition/erosion formulation is applied. Sediment erosion begins after the current reaches a steady state, during which eddy viscosity is of a parabolic shape. Simulations are then performed for a three-hour duration, which gives a significant change in depth due to erosion in the case of the strongest wave.

Figure 2 presents results of different modeling approaches for the bottom sink/source terms in different wave conditions. In the case of ${\tau }_w=0.85$ m s⁻¹ (Figure 2a), the wave-induced shear stress is relatively weak and different modeling approaches give almost the same results. As ${\tau }_w$ increases, the difference among cases using different approaches becomes visible, as shown in Figure 2b. In the case when the wave-induced stress is the largest, as shown in Figure 2c, a significant overestimate of SSC can be found using the original erosion/deposition formula. Figure 2c also shows that the modified erosion/deposition formula greatly eliminates the overestimates, which demonstrates the importance of applying the modified erosion/deposition formula when suspended sediment modeling is conducted on the fixed-grid framework. Only little differences are found between cases of the original erosion/deposition formula on moving grid and of the modified erosion/deposition formula, which is due to the difference of the bed shear stresses induced by changing depth in the moving-grid case. Comparison of SSC profiles in different wave conditions in Figure 2 shows that SSC profiles dramatically increase as ${\tau }_w$ increases. This is due to the exponential formula (Equation (21)) employed for bed erosion, which is potentially ill posed when applied to the strong wave case on fixed grids without the present modified erosion/deposition formula.

A further run with the strongest ${\tau }_w$ is then conducted so that simulations reach the steady state. Although it is not a realistic situation that the critical shear stress remains the same when the bed is being eroded (i.e., the depth-varying ${\tau }_{CE}$ should be considered due to mud consolidation), the objective is to examine the potential difference among different modeling approaches in a long-time run using a simple model setting. Figure 3 presents the time series of SSC in the bottommost cell from time = 0 to 24 h, which is the time when all simulated SSCs almost reach the steady state. It can be seen that both the moving-grid and the modified erosion/deposition cases nearly reaches a steady state at time = 10 h, while the case of the original erosion/deposition formulation with fixed grids reach a steady state at time = 24 h. The growth rate of SSC in the case of the original erosion/deposition formulation is faster than the other two cases throughout the simulation, and without any connection to the properties of the bed material in the original erosion/deposition formulation, SSC in the fixed grid simulation exceeds the prescribed dry density (530 g L⁻¹ here) at approximately 10 h.

In this section, we show that, for a simple 1DV model, slight changes in the wave condition may dramatically increase the quantity of eroded sediment. In such a case, implementation of the original erosion/deposition formula disregarding the change in depth can lead to unrealistically high SSC in the modeling result due to the ill-posedness when applying the exponential law to model erosion. Moreover, we demonstrate that, on the non-moving grid, this issue can be fixed by applying the present modified formula, which is able to give SSC results similar to those with moving grids. Therefore, when suspended sediment transport modeling is incorporated with the fixed-grid circulation model, it is recommended to apply the present modified erosion/deposition formula to obtain the stabilized SSC. In the next section, an example of field modeling is presented.

3.2. Three-Dimensional Sediment Transport Modeling of San Francisco Bay, California, USA

In the second example, we present three-dimensional (3D) modeling of hydrodynamics and sediment transport in San Francisco Bay using the Stanford unstructured coastal hydrodynamics model, SUNTANS. In San Francisco Bay, a distinct channel extends throughout the bay along its thalweg and incises shallow shoals (see Figure 4a). This results in the intriguing dynamics of sediment transport subject to combined wave and tidal forcing, which has been extensively studied in past few decades (e.g., [28,38,39,40,41]). Here, we focus on the comparison of modeled SSC between the conventional and present modified formulation at North San Francisco Bay (North Bay) (i.e., San Pablo Bay and Suisun Bay, see Figure 5) using the validated SUNTANS model, with an emphasis on the consequence of wave forcing. The domain setup of Chua and Fringer [42], who studied tidally driven salinity in San Francisco Bay during January of 2005, is employed here. The same setup was also employed to study hydrodynamics due to waves and tidal currents at South San Francisco Bay during September of 2009 by Chou et al. [43] and the sediment transport by Chou et al. [5]. All these studies showed good agreement compared to field observations. More details of the model and simulation setup can be found in these studies. As we focus on the effect of the modified formulation for modeling erosion and deposition, in what follows, only setup conditions that are relevant to the present results are addressed.

The computational domain spans between the Pacific Ocean and the Sacramento-San Joaquin River Delta, as shown in Figure 5. Open boundaries are at the Pacific Ocean, Sacramento River, and San Joaquin River. The average grid resolution based on triangular cell lengths is 50 m. In the vertical, the z leveled grid has a maximum of 60 layers in the deepest portion of the domain. The vertical resolution is refined near the surface with a grid stretching ratio of 10%, resulting in a minimum vertical grid spacing of 0.29 m. The total number of cells in the horizontal is approximately 80,000. The flow is forced by specifying free-surface elevation along the Pacific Ocean boundary using the data at Point Reyes (see Figure 5), which is obtained from the National Oceanic and Atmospheric Administration (NOAA) Center of Operational Oceanographic Products. Flow rates at the Sacramento and San Joaquin Rivers are given by freshwater inflow estimates from the DAYFLOW program of the California Department of Water Resources. A time step of 10 s is employed to ensure stable horizontal advection of scalars. Waves are generated due to winds, and their transport is modeled using the wave-action approach following Booij et al. [44] on the same grid with the same time step as the hydrodynamics model. Wind speed and direction are reconstructed as each grid cell by interpolating wind data from five NOAA wind stations, as shown in Figure 5. More details of the wave model and simulation setup can be found in Chou et al. [43]. For the sediment transport model, according to the field study by Manning and Schoellhamer [45], two size classes are used to model a slowly- and fast-settling components, respectively. For the bed, a multi-layer bed sediment model [5] is used. In this model, the bed is vertically divided into three layers, and each layer is assigned with different erodibility ($E_{b,0}$), critical shear stress for erosion (${\tau }_{CE}$) (see Equation (18)), height $h_L$, and dry density $\rho bd$. Moreover, consolidation rates $T_c$ are assigned to act as a mass exchange between bed layers. More details of the bed model can be found in Chou et al. [5]. Properties of the bed sediment are summarized in Table 1. Resuspension of the bed sediment into the water column is modeled in the way that once the bottom stress exceeds the critical shear stress for erosion, quantities of resuspended sediment are calculated by Equation (21) and divided into each size class according to prescribed fractions. The fraction of each size class, as denoted by $R_j$ (j in the class index), is empirical and needs to be calibrated with the field measurement. According to Chou et al. [5], with $w_{s,1}=0.1$ (fine) and $w_{s,2}=2$ mm s⁻¹ (coarse), $R_1$ = 5% ($R_2$ = 95%) in regions where the depth < 5 m and $R_1$ = 30% ($R_2$ = 70%) in regions where the depth ≥ 5 m give the best agreement with the field data. Here, following previous studies, we use these setup conditions with the original and modified resupension formulas for comparison.

Figure 6 presents the time series of the root mean squared (RMS) wave height, $h_{w,rms}$, measured near the Port Chicago wind station, as indicated by the cross in Figure 4b. During the summer, winds coming from the Pacific Ocean result in persistent wave events in the Bay. As shown in Figure 6, the wave event persists during yeardays 255–256. Here, we choose the time point of 4:30 p.m. on the yearday 255 (13 September 2009) for the following discussion. This time point is chosen also because it is the time when SSC predicted by the original resuspension formula starts to deviate from those by the modified formula and becomes unrealistic. The corresponding surface elevation and the resulting distribution of the RMS wave height along with the wind field are shown in Figure 4b,c, respectively. During the time, as shown in Figure 4c, modeled waves are relatively strong at the Suisun Bay without significant breaking, which results in strong near-bottom orbital velocities. As shown by the transect in Figure 7a, particularly in regions of low water depth, this leads to high near-bottom SSC, which can further transport into the deep-water channel [5]. However, without modification, the exponential relationship given by Equation (21) can easily lead to extremely high sediment concentration near the bottom, which subsequently makes the SSC in the water column unrealistically high, as shown in Figure 7b. A comparison between Figure 7a,b clearly shows that the modeled near-bed SSC can be stabilized by the present modification. The temporal evolution of the horizontal distribution of the near-surface SSC in the case with the original resuspension formulation is presented in Figure 8. It can be seen that the modeled SSC begins to grow significantly and forms certain spots of the high concentration (see Figure 8a,b). The subsequent dramatic growth during the persistent wave events results in unrealistically high SSC in the channel (see Figure 8c,d), which can be eliminated by applying the present modification for the suspension formula. Figure 9 shows comparison of time histories of the point measurement for the near-surface SSC between cases with modified and original resuspension formulas, along with the case without the wave effect. First, compared with the case without waves, Figure 9 shows that considering the wave effect (with the modified formula) significantly enhances SSC in the water column, indicating the dominance of waves in resuspending bed sediment. In the present case, the strongest wave effect in the case with the modified formulation almost doubles the SSC compared with the case without waves (during yeardays 256–257). Moreover, without any stabilizing effect for the resuspension formula, Figure 9 also shows an unrealistic growth in SSC after the yearday 256, which demonstrates the importance of the present modification when modeling the wave-induced resuspension.

In fact, it should be noted that a more detailed field survey for the realistic bed properties, such as the erodibility as a function of the erosion depth, may also prevent the unstable growth of the near-bed SSC. For example, supply of fine sediment, such as silt and clay, is always limited, and a realistic bed parameter should be able to immediately limit the bed supply when the stress is too high. However, a thorough field campaign to obtain proper parameters is expensive. In addition, complete and representative datasets are usually unavailable. Here, we present a simple modification that considers the realistic SSC into the water column during resuspension and automatically stabilizes the SSC modeling in the estuarine circulation model. This, as demonstrated here, is of critical importance when waves are considered.

4. Conclusions

The scalar transport equation in conjunction with the hydrodynamics model has emerged as a standard approach to modeling the transport of suspended sediment. As a result, the near-bed source/sink terms attributed to erosion and deposition play instrumental roles in model capability. Our results demonstrate that, regardless of the type of formulas used to describe erosion and deposition, current implementations that disregard volume conservation could lead to an overestimation of SSC. Thus, we derived formulas modified with respect to moving control volume, i.e., changes in water depth due to sediment erosion and deposition. The modified formulation constrains the near-bed SSC such that it is always less than the prescribed dry density. In addition to giving more realistic SSC, the modified formula guarantees model stability when the stratification effect is taken into account. Numerical examples are presented using the 1DV suspension model, the results of which show that the original erosion/deposition formulation can give unrealistically excessive suspension on fixed grids. However, this can largely be overcome using the proposed modification. As in the present model set-up, the proposed concept of modified erosion/deposition formulation is particularly important in cases of strong bed stress. A typical example is a shallow water environment in which surface gravity waves can induce a strong shear on the bed.

In the first test case, using a 1DV model, we simply demonstrate a significant deviation of the conventional fixed-bed modeling result from the realistic moving-bed result when modeling the resuspension in a turbulent water column subject to wave forcing. In particular, the former case can give unrealistic results in that the near-bed SSC exceeds the prescribed dry density of the bed material. We also show that this issue can be fixed by the present modification for the erosion formula. In the following test case of 3D field modeling of San Francisco Bay, we show that the present modified formulation can significantly suppress the resuspension of the bed sediment. This is of particular importance in shallow estuaries, in which the wave effect plays an important role in resuspending the bottom mud. When incorporated with the power- or exponential-law to model sediment resuspension, the addition of the wave-induced shear stress can easily lead to unrealistically excessive erosion. In this case, the present example demonstrates an effective stabilization for the modeled SSC using the proposed modification for the erosion/deposition formula.

Author Contributions

Conceptualization, Y.-J.C.; Methodology, Y.-J.C., Y.-C.S., Y.-H.S. and C.-J.C.; Software, Y.-J.C. and Y.-C.S.; validation, Y.-J.C., Y.-C.S. and Y.-H.S.; Formal analysis, Y.-J.C.; Investigation, Y.-J.C.; Data curation, Y.-H.S. and C.-J.C.; Writing–original draft preparation, Y.-J.C.; Writing–review and editing, Y.-J.C.; Visualization, Y.-H.S.; Supervision, Y.-J.C.; project administration, Y.-J.C.; funding acquisition, Y.-J.C.

Funding

This research was funded by Ministry of Science and Technology (MOST), Taiwan, Civil and Hydraulic Engineering Program under Grant No. 106-2628-E-002-016-MY4, the Ministry of Education (MOE) Higher Education Sprout Project in Taiwan (NTU Core Consortium: 108L880206), and National Taiwan University career development project (NTU-CDP-108L7720). The APC was funded by Ministry of Science and Technology (MOST).

Acknowledgments

The computation was performed on the Dell PowerEdge M620 supercomputer cluster located at the Center for Advanced Study in Theoretical Sciences, which is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Figures and Table

Figure 1. Schematic plot to show modeled deposition (a); real deposition (b); modeled erosion (c); and real erosion (d).

Figure 2. Vertical profiles of SSC for τw=0.85 (a), 1.32 (b), and 1.74 N m−2 (c) in the cases of modified erosion/deposition formulation with fixed grids (black solid line), original formulation with moving grids (black dashed line), and original formulation with fixed grids (gray solid line) in different wave conditions.

Figure 3. Time series of SSC at the bottommost cell in the cases of modified erosion/deposition formulation with fixed grids (black solid line), original formulation with moving grids (black dashed line), and original formulation with fixed grids (gray solid line).

Figure 4. Model bathymetry (in m below the mean higher high water (MHHW)) of North San Francisco Bay (a), following Chua and Fringer [42], in which the white contours indicate a depth 5 m to separate shallow-water shoals (depth < 5 m) from the deep channel (depth > 5 m), and modeled surface elevation (b) and RMS wave height (hw,rms) (c) at 4:30 p.m. on 13 September 2009. The cross in (b) indicates the location of the point measurement of SSC and wave height. Distances are relative to the point measurement marked by the cross in (b). The gray solid line (AB¯) in (b) indicates the transect of SSC profiles presented in Figure 7.

Figure 5. San Francisco Bay coastlines along with the locations of wind stations (black triangles) and the tidal gauge (black square).

Figure 6. Time series of the modeled root mean squared (RMS) wave height (hw,rms) measured at the point indicated by the cross in Figure 4b). The gray dashed line indicates the time point of 4:30 p.m. on 13 September 2009, which corresponds to the time point in Figure 4, Figure 7, and Figure 8a.

Figure 7. Modeled SSC along the transect shown in Figure 4b at 4:30 p.m. on 13 September 2009 in cases with the modified resuspension formula (a) and with the original resuspension formula (b).

Figure 8. Representative snapshots of the near-surface SSC at 04:30 p.m. (a), 05:30 p.m. (b), 06:30 p.m. (c), and 07:30 p.m. (d) on 13 September 2009 predicted using the original resuspension formula without the present modification.

Figure 9. Comparison of time series of the modeled near-surface SSC measured at the point indicated by the cross in Figure 4b) using the original resuspension formula with (gray solid) and without (black dashed) waves and the modified resuspension formula with waves (black solid).

Table 1

Parameters of the multilayer bed model. Here, ${\rho }_{bd}$ is the dry density, $T_c$ is the consolidation rate, $h_L$ is the layer height, and the other parameters are defined in Equation (21).

Layer No.	${\mathit{\rho }}_{\mathit{bd}}$ (g m−3)	${\mathit{\tau }}_{\mathit{CE}}$ (N m−2)	$\mathit{\alpha }$	${\mathit{E}}_{\mathit{b}\mathbf{,}\mathbf{0}}$ (g m−2 s−1)	$\mathit{\beta }$	${\mathit{T}}_{\mathit{c}}$ (g m−2 s−1)	${\mathit{h}}_{\mathit{L}}$ (m)
1	75,000	0.1	4.5	0.01	1	0.0002	0.1
2	530,000	0.4	4.5	0.01	1	0.0002	0.5
3	1,200,000	1.2	4.5	0.01	1	0.0002	4.0