Using a Holistic Modeling Approach to Simulate Mud‐Induced Periodic Stratification in Hyper‐Turbid Estuaries

This study focuses on a holistic modeling approach, in which water, fluid mud, and immobile mud are all calculated by only one set of equations. To integrate the immobile mud into this concept, a holistic transport equation including sediment transport and consolidation is developed. In some estuaries, extensive deepening and dredging resulted in tidal deformation and sediment import to such extent, that hyper‐turbid conditions developed. Recent measurements from the Ems estuary show that the locations of interfaces between water, fluid mud, and consolidated mud vary during a tidal cycle. Conditions are varying from fully mixed to stably stratified. As a suitable case study for the holistic model, a 1D vertical numerical simulation of the Ems has been set up, which is able to qualitatively reproduce the observed vertical velocity, concentration, and velocity shear profile. The simulation shows mud‐induced periodic stratification.


of 12
During ebb tide (position C) the flow velocities are usually not as large as during flood tide due to tidal asymmetry. Accordingly, the turbulent mixing is weaker compared to flood tide such that it cannot destroy the layered structure entirely. Figure 1c shows the sediment concentration in the Ems estuary from optical backscatter measurements of Becker et al. (2018a). The sensor casts were collected every 30 min and the resulting vertical profiles are interpolated over time.
In order to investigate the impact of future measures to improve the situation in hyper-turbid estuaries, numerical models should be able to simulate the actual behavior of the system. When this is possible, the impact of planned measures can be studied using a numerical model before they are implemented.
Up to now, the interaction between fluid mud and tidal dynamics is seldom considered in hydrodynamic-numerical models for large-scale applications. The simulation of fluid mud dynamics is complex because of specific rheological properties of mud and a wide range of interactions between sediment characteristics and flow behavior. Le Normant et al. (1993) used the 3D free surface flow model TELEMAC-3D and included a sediment module to simulate cohesive sediment problems. The fluid mud dynamics were modeled by a depth-averaged approach, where fluid mud was represented by a constant, high viscosity. Le Normant (2000) applied this model for the Loire Estuary with a combination of a fluid mud module, erosion, deposition, and a vertical 1D bed consolidation module. Another fluid mud module is available in Delft3D, which calculates a separate mud layer with a constant sediment concentration . However, those models neglect the rheological shear-thinning behavior of mud as a function of velocity shear and sediment concentration. Wehr and Malcherek (2012) simulated the complex non-Newtonian behavior of fluid mud with an isopycnal numerical model. This model concept can reproduce stably stratified flow conditions but not vertical mixing processes. Le Hir et al. (2001) introduced a continuous modeling approach, which combined water and mobile mud. Roland et al. (2012) developed a fluid mud module based on that continuous approach, which was coupled with the 3D flow models SHYFEM and SELFE. The model is able to reproduce stable stratification as well as vertical mixing of fluid mud. However, the immobile mud cannot be modeled and a detailed application and validation for an estuary are still needed. There is actually no 3D numerical model available which is able to simulate the described behavior of the sediment concentration 3 of 12 during a tidal cycle in the Ems. Therefore, as a first step, 1D vertical (1DV) models should be developed in order to understand the fundamental processes involved in that complex system.
In this study, a holistic approach is presented which uses the same set of fundamental balance equations for the water column as well as for the mud layer. The model can be regarded as an extension of the continuous model approach as it goes one step further by including the immobile mud. Deposition, erosion, and consolidation are inherent parts of the holistic equations.

Materials and Methods
The numerical model consists of four differential equations: The momentum balance, the sediment transport equation, and the two-equation k- E turbulence model. The four equations are implemented using the finite volume method.
Assuming that the main processes take place in the vertical direction only, a 1DV numerical simulation has been set up. This means that horizontal divergence of transports is not considered. Furthermore, earth rotation, as well as internal pressure gradients due to horizontal density gradients are neglected. The suspended sediment is not treated as a separate phase but represented by a volumetric concentration and a settling velocity. The simulation does not include salinity, as salinity gradients in the Ems estuary are much smaller than those of sediment concentrations.

The 1D momentum equation in vertical (z) direction is
. Combining rheological and turbulent viscosity was introduced by Le Hir et al. (2001). The rheological viscosity can be derived from measurements of fluid mud with a rheometer. The use of a rheological model based on specific measurements was implemented in a numerical simulation for the first time by Knoch and Malcherek (2011). The rheological viscosity is large when sediment concentrations are large and velocity shear is low. The turbulent viscosity is calculated by the turbulence model. In estuaries, a high turbulent viscosity results from strong turbulence due to vertical shear. By contrast, a high rheological viscosity suppresses turbulence and can lead to laminar flow conditions. The adapted 1D sediment transport equation reads The yield stress increases with concentration, due to the increase in the number of interactions and the amplification of their intensity (Coussot & Piau, 1994 ) has a non-zero Bingham yield stress, which is not correct. For this reason, we adopt the power relation where in the case of   0 s E , there is no yield stress. There exist many yield stress measurements of mud samples for different densities and locations. Figure 2a compares the yield stress from different estuaries as a function of the volume fraction of solid particles. A least squares fit was applied to the data shown. We set  4 E r and achieved   5 The dependence of the Bingham viscosity on the solid content is The rheological viscosity increases with higher volume fractions and decreases with increasing shear rate. The latter occurs due to the shear-thinning behavior of fluid mud. The rheological viscosity as defined by Equation 3 is a possible parametrization for mud in the Ems estuary and is valid for    0 1 s E .

Turbulence Model
Turbulence models are used to calculate the turbulent viscosity in a fluid that affects several properties, for example, the velocity in the flow. Generally, in a stationary shear flow without sediment, a logarithmic velocity profile develops. Very close to the wall, in the viscous sublayer, the velocity profile is however linear.
for high concentrations. The turbulence model by Wilcox (1988) and adapted from Chmiel et al. (2020) is

A Holistic Sediment Transport Equation
Sediment transport and consolidation affect each other. If sediment particles settle due to a decrease in turbulence, the sediments form a highly concentrated suspension near the bottom. If shear forces remain low, the suspension starts to consolidate. If shear forces are high, sediment that is already consolidated can be liquefied and entrained into the water column. Consequently, to account for the dependence of both processes, they will be modeled together in one single differential equation. To do so, the Gibson consolidation equation needs to be rearranged such that it has the same mathematical shape as the classical sediment transport equation. This means that it consists of an advective and a diffusive term for the sediment concentration E c . The one-dimensional consolidation equation from Gibson et al. (1967) is This new equation now has the shape of the classical sediment transport equation and, as a result, leads to a formulation for a consolidation settling velocity and a consolidation diffusivity. The consolidation settling velocity is s c c f s s w w k The consolidation diffusivity is The permeability f E k is derived from Kozeny (1927) and Carman (1937)

From Settling to Consolidation
The next step is to establish a relation between settling in the water column and settling due to consolidation. Fluid mud can be described as fine-grained sediment particles in a water suspension. As sediment has a higher density than water, the particles settle within the water column. For a single, spherical particle with the diameter E d , Stokes' formula for the settling velocity is valid. However, falling particles generate a return flow and influence other particles in the vicinity. Thus, the settling velocity of the suspension decreases by a factor   (1 ) s E . When comparing Equations 8 and 10 similarities become obvious. In order to achieve a transition for both settling velocities, a blending is suggested such that which is typical for high concentration suspensions. Winterwerp (2002) states three reasons for hindered settling. First, the return flow of a particle, which reduces the settling velocity of surrounding particles. Second, the increase of the viscosity of the suspension and, third, enhanced buoyancy due to increased concentration. All effects reduce the settling velocity with increasing sediment concentration, which leads to 7 of 12 . 1 2.5  As illustrated in Figure 3b,  s E decreases with increasing velocity shear and decreasing turbulent viscosity. Typically,  s E is a function that is derived from second-order turbulence models, depending on shear and stratification. Here, the turbulent Schmidt number is also a function but not derived from second-order closure but from adding the erosion diffusivity to the turbulent viscosity. Consequently, the relation between mixing of momentum and mixing of masses is modified. Further research is needed to verify this assumption. This modification also influences the buoyancy term, which is typically not only present in Equation 4 but also in Equation 5. When the buoyancy term is added, another parameter that is dependent on the Schmidt number function arises. Again, as the relation between mixing of momentum and mixing of sediment was modified, this is also subject to further research.

Simulation of the Ems Model
The aim of the simulation is to reproduce the processes occurring in the Ems estuary, especially to simulate mud-induced periodic stratification and to show that the model includes the immobile bed. The 1DV model solves the differential equations for velocity, sediment concentration, and the turbulent quantities E k and  E . The model Equations 1, 2, 4, and 5 are discretized using the finite volume method. The water level E h and the gravitational acceleration E gJ (where  dh E J dx ) are predefined. In order to account for the time-dependent water level, tidal constituents M2 and M4 are considered, as they are the most dominant at the study site.
The M4 component has a phase shift that yields the characteristic asymmetrical profile of the Ems. The flood phase is shorter and exhibits higher velocities compared to the ebb phases. The vertical resolution is   0.237 m E z ; the grid is equidistant except for the last cell, which is adapted to the current water level. We use Equation 3 for the rheological viscosity, Equation 11 for the settling velocity, and Equation 13 for the holistic diffusivity. An erosion mixing length of    3 6.7 10 m c E l yields the best model fit to the observations from the Ems estuary (Table 1).
The 1DV model calculates the velocity, the sediment concentration, and the turbulent parameters over the vertical simulation domain. The total simulation time is 5 days, the time step  E t is 10 s. The simulation reaches a steady state after an initial time period of three tidal cycles. The initial concentration in each cell is . This value is chosen to be small, to reduce the flow velocity at the free surface. The model has been implemented in MATLAB. The simulation results are most sensitive to the vertical resolution, the settling velocity, the rheological viscosity, the erosion mixing length, and the initial sediment concentration. Figure 4a shows the simulation results of the vertical concentration profile over a time period of 19 hr in analogy to Figure 1c. Figure 4b compares the concentration, velocity, and velocity shear profiles for the positions A, B, and C for measurements and simulation. The water level periodically rises and falls similarly to the actual water level. At position A the sediment is distributed over the entire water column and the immobile mud layer is the thinnest among all three positions because of enhanced erosion during flood tide. The velocities at position A are high and directed upstream, and the highest shear is located close to the immobile mud. The simulation deviates from the linear concentration profile from the measurements. A possible explanation is that in the simulation only M2 and M4 tidal components are considered. As a consequence, the beginning of flood tide sets in less abruptly and leads to less acceleration and less mixing Potential of a dissipation rate 1 E s  ) are slightly higher than the actual values. The velocity profile at the same position consists of lower velocities in the fluid mud layer and higher velocities above the lutocline. The velocities in the fluid mud are overestimated by the simulation, which could be due to too low rheological viscosities. The velocity shear yields similarly high values above the upper and the lower lutocline. The field observations show higher velocity shear above the upper lutocline but the general behavior is well captured. Accordingly, the simulation is able to reproduce the stratification of fluid mud. During ebb tide (position C), the velocities are directed downstream and the position of highest shear is right above the lower lutocline. The measurements deviate from that, as the upper lutocline still exists at position C. The measurement data does not provide insight into velocity shear data close to the bottom, which is probably lower than the calculated values. The too rapid destruction of the lutocline with the onset of ebb tide is probably due to several effects. First, the ebb tide starts too early compared to the measured data, causing the flow velocities to be too high. Second, the neglect of flocculation processes could also play a role: Macroflocs, which have higher settling velocities, could form during slack water and stabilize the lutocline at the beginning of the ebb phase. The simulation results should also be considered from the point of view of the 1D model: A vertical 1D model is only capable to simulate vertical transport processes (settling, entrainment, deposition, and erosion). It is not valid when horizontal transport processes are dominant in an estuarine system.

Results and Discussion
Although there are more advanced rheological models, turbulence models, or settling velocity approaches than those presented here, the simulation is able to qualitatively reproduce the mud dynamics in the Ems estuary. These parameters, which are usually determined by field or laboratory measurements, are subject to large uncertainties. The 1D model can therefore be considered as a useful tool to validate these parameters, which can then be implemented in 2D and 3D simulations. The total number of parameters used in our simulation was kept to a minimum, and the simulation of the Ems presented here does not require a separate fluid mud or sediment module. Overall, the simulation shows mud-induced periodic stratification and the results are in fair agreement with the actual data.

Summary and Conclusions
A model approach has been presented which consists of four differential equations to calculate the flow velocity, concentration, and turbulence profile in a hyper-turbid estuary. The Gibson equation for consolidation was rearranged to the same structure as the sediment transport equation. This procedure allows to include the simulation of sediment transport processes from settling to consolidation in one single equation. The model is valid in the water column, the moving mud as well as in the immobile bed. In that way, the model can be regarded as holistic. In classical models, the described layers are separated, and connected with empirical parametrizations.
A simulation of the hyper-turbid Ems estuary is presented and the simulation results show that the model can reproduce mud-induced periodic stratification. The simulation results are comparable to recent measurements in the Ems estuary for the concentration profile, flow velocity, and velocity shear. The latter is an indicator for the production of turbulence and, on that account, turbulence itself. All quantities are also calculated in the immobile bed, where the flow velocity and turbulence are zero and concentrations are very high.
The holistic modeling approach is applicable to various problems, for example to studies of sediment settling. To simulate a settling column for fine sediments, the free surface gradient in the momentum balance has to be set to zero. A homogeneous concentration can be defined at the initial time step. With time, the sediment settles and consolidates. This means that the model setup can be used to calibrate model parameters for sediment settling and consolidation.

Data Availability Statement
V1.1 of the model code of the 1DV model used for the simulation is preserved at Malcherek and Schmidt (2021). The measurement data used for comparison is available through Becker et al. (2018b).