Wave-Created Mud Suspensions : A Theoretical Study

We studied wave-created high-density mud suspensions (fluid mud) using a one-dimensional water column (1DV) model that includes k-ε turbulence closure at a high vertical resolution with a vertical grid spacing of 1 mm. The k-ε turbulence model includes two sediment-related dissipation terms associated with vertical density stratification and viscous drag of flows around sediment particles. To this end, the calibrated model reproduces the key characteristics (maximum concentration and thickness) of fluid mud layers created in laboratory experiments over a large range of wave velocities from 10 to 55 cm/s. The findings demonstrate that the equilibrium near-bed mud concentration, Cb, is solely determined from the balance between erosion and deposition fluxes, whereas the thickness of the fluid mud layer is mainly controlled by sediment-induced density stratification, which dissipates turbulence and hence eliminates turbulent sediment diffusivity at the top of the fluid mud layer, the lutocline. Our model stands in contrast to those that suggest that upward sediment diffusion is close to zero at the interface between the fluid mud layer and the overlying fluid. Instead, our model suggests that the upward diffusive flux of fluid mud flows peak at the lutocline and is compensated for enhanced settling fluxes just above it. Our model findings also support the existence of the gelling-ignition process, which is critical for the development of fluid mud beds in modern sedimentary environments.


Introduction
Fluid mud is the mixture of clay, silt and very fine sand particles, organics, and water at concentrations close to the tightest possible packing of mud floccules, termed the gelling point, C gel (Figure 1).Fluid mud layers are often only a few centimeters thick [1,2], but they can reach much greater thicknesses of up to several meters [3][4][5].
The behavior of suspended mud in waterways has been the focus of intense research, in part, because mud is typically enriched in toxins (e.g., heavy metals) and because suspended mud may hinder the marine navigation of vessels.Improved predictability of mud movement is critical for coastal, estuary and river management, and for better understanding of the dynamics and biogeochemistry of fluvially derived mud on continental shelves, which plays an important role in the nutrient and carbon cycles in the oceans [6].
Suspended mud may accumulate gradually by sedimentation from the water column, either as individual clay particles or as floccules.In most settings (not adjacent to deltas), this is a slow process, as the sediment concentration in the water column is usually small (<1 g/L = 1 kg/m 3 ), and the level of flocculation and settling velocity of individual mud particles is relatively low.Another transport mechanism is the movement of concentrated, near-bottom mud layers, often referred to as fluid mud, under the influence of waves and gravity [7].Bed shear stresses induced by short-period waves such as wind waves and swell can be several orders of magnitude higher than for other slowly varying flows such as tides [8].Hence, wave-seabed interaction is a prime mechanism creating fluid mud, and may result in high mud transport rates during a short period of time, as sediment concentrations in fluid mud may be up to ~100 g/L.This transport mechanism is deemed responsible for the rapidly deposited mud accumulations often observed after storms [9].Fluid mud layers of several meters in thickness have been documented on the continental shelf [3] and estuaries [4].Such thick fluid mud layers are often attributed to the existence of an external source (e.g., river discharge), a vertical sediment-flux convergence (sedimentation rate < consolidation rate), strong-tidally induced turbulence, and a horizontal flow convergence [3].When the concentration of floccules becomes high enough, settling floccules start to hinder each other in their movements, a process known as "hindered settling" [11].Our previous work [12] analyzed the bottom-boundary equation governing the flux of cohesive mud across the seabed.Findings suggest that higher bed shear stresses can trigger a situation in which a positive feedback between hindered settling and the net erosion flux leads to rapidly increasing suspended-mud concentrations towards the gelling point, a mechanism that we termed "gelling ignition" (Figure 2).A recent three-dimensional model application [13] confirms and discusses the existence of this mechanism.The viscosity of fluid mud substantially increases at higher concentrations near the gelling point [14] and fluid mud becomes governed by non-Newtonian effects [15].In addition, the settling flux increases due to permeability effects and induces self-weight consolidation [16,17].
Reported gelling concentrations (i.e., suspended sediment-concentrations at the gelling point) are in a range of 40-120 g/L [11], noting that ambient turbulence can shift the gelling point to substantially higher concentrations [18].These gelling concentrations correspond to mud-seawater bulk densities in a range of 1050 to 1100 kg/m 3 .Owing to weight excess, such fluid mud layers can facilitate density-driven flow [11], transporting suspended mud across the seabed.An interesting aspect of fluid mud is its ability to move along vast distances on extremely mild bottom slopes [19].For instance, slurries travel distances over 40 km at a bottom gradient of 1:10,000 in the Sanmenxia Reservoir in China [20].
Early laboratory experiments of mud erosion by waves [2] considered forcing periods in a range of 1-2 s, using kaolinite and estuarial-mud beds with consolidation times in a range of 2 to 15 days.Wave-induced bed shear stresses were relatively low (0.11-0.2 Pa) and created thin (0.5-1 cm) fluid mud layers of relatively constant concentration (~100 g/L).These authors [2] suggested that the erosion flux for cohesive sediment under waves can also be expressed by the classical Ariathurai-Partheniades formula, which has been derived for steady flows, with use of a small value of the bed shear strength (bed shear resistance to erosion) of only 0.02-0.1 Pa.Wave-induced bed liquefaction facilitates this lowering of bed shear strength [21].Only few systematic laboratory experiments considered wave-created mud suspensions [22,23].In these experiments, suspended fluid mud layers were created by oscillatory flows of periods between 3 and 8 s and maximum velocities in a range from 15 to 60 cm/s over a sediment bed (silica flour) with a mean particle size of D = 20 µm.The predominantly silt-sized mixture contained smaller fractions of ~10% clay (D < 3.9 µm) and ~20% sand (D > 63 µm).Findings of this laboratory study serve as a benchmark for the numerical model predictions presented here.
Fluid mud of concentrations of up to 150 g/L was observed on the Po prodelta off the Italian coast of the Adriatic Sea [1].Its formation coincided with storms that created wave velocities of 20-50 cm/s and bed shear stresses of 1-3 Pa.The observed high-density mud suspensions had thicknesses of 5-8 cm.Other studies [24,25] report the observation of liquefaction of the seabed prior to wave-induced erosion and the subsequent formation of fluid mud layers (~20 cm thick).
The modeling of fluid mud is challenging for many reasons.The dynamics of fluid mud transport involve a variety of complex physical mechanisms, including for example, gravity-driven transport, turbulence modulation, flocculation, hindered settling, non-Newtonian rheological behavior [16], and sediment fluxes across the base of the fluid mud layer being modified by both consolidation and wave-induced fluidization (see Figure 1).Typical coastal models, such as DELFT3D [26] or NOPP-CSTM [27] are not designed to resolve the vertical structure and dynamics of high-density mud suspensions.Hence, processes that occur within this layer need to be parameterized to derive quantities such as bottom drag coefficient and sediment transport rate.
There are few numerical-modelling frameworks for cohesive sediment transport that vertically resolve thin fluid mud layers.Whereas earlier one-dimensional water-column (1DV) models for mud suspensions [28,29] incorporated the damping of flow turbulence due to sediment-induced stratification, which is a crucial mechanism for the formation of fluid mud, the rheology of sediment was neglected.Following advancements in the understanding of multiphase noncohesive sediment transport [30,31] and physical-based description of both floccule dynamics [32,33] and consolidation [17], the first comprehensive numerical models for cohesive sediment were only developed recently [34,35].Despite this, only a few of these model experiments explicitly addressed the formation of fluid mud with high suspended sediment concentrations of ~100 g/L.
This work focuses on the creation of fluid mud layer via wave-induced erosion of bed sediment.To this end, we present an updated 1DV model and discuss parameter configurations that give results similar to laboratory experiments [22,23].Other processes by which fluid mud flows form, such as bulk erosion/instabilities of the consolidated mud bed via liquefaction and subsequent fluidization, are not considered in this work.

Materials and Methods
This section describes the physical equations governing the dynamics of mud suspensions, algorithms of erosion-deposition fluxes, and research approaches used in the study.The latter consist of a stand-alone analysis of sediment fluxes across the seafloor and applications of a one-dimensional water-column model, known as 1DV model.We justify the use of a one-dimensional model by the fact that surface gravity waves impacting the seafloor are classified as shallow-water waves (also known as long waves) that operate on a horizontal wavelength much longer than the total water depth.Hence, relatively large areas of the seabed are expected to respond in unison to the wave forcing, and there is little horizontal variability on such scale.The presence of shallow-water waves also implies that wave-induced vertical velocities are extremely small and negligible in the dynamics of fluid mud.

1DV Model: Wave Forcing and Momentum Equation
The water column model, described in the following, is applied at a vertical grid spacing of ∆z = 1 mm extending a total distance of 50 cm from the seafloor.The model is forced by an oscillatory ambient horizontal flow of the form of: where u a is the maximum speed of the oscillatory horizontal flow, t is time, and T is wave period, set to 6 s.Variations of T did not significantly alter the results.The total simulation time of experiments is 30 min.Within the wave boundary layer, which has typical thicknesses of 3-10 cm [8], the horizontal flow changes under the influence of bottom friction and vertical momentum diffusion according to: where u is the deviation of the horizontal velocity from the prescribed barotropic flow (1), z is the vertical coordinate, A z is total vertical viscosity, and ϕ = C/ρ s is sediment volume concentration with C being suspended sediment concentration and ρ s = 2650 kg/m 3 (sediment density approximated as quartz).From (2) it should be clear that, apart from temporal variability, the simulated wave-induced horizontal velocity field, given by u w + u , strongly varies vertically due to diffusive and frictional influences.Note that (2) assumes a plane seafloor; that is, any slope effects are not considered in this study.Total vertical viscosity in (2) is given by: where ν m = 1.4 × 10 −6 m 2 /s is kinematic molecular viscosity, ν t is turbulent viscosity, derived from a k-ε turbulence closure, detailed below, and ν e denotes non-Newtonian rheological effects.The latter is calculated from a Bingham-like formulation [36] according to: where ν ∞ , ν 1 , and γ are parameters, and ∂u/∂z is the shear strain rate.The time scale γ determines the magnitude of shear-thinning effects [37].This study adopts the parameter setting (ν ∞ = 0.003 m 2 /s, ν 1 = 3 m 2 /s and γ = 1000 s) of previous studies [35,36].Note that in the experiments discussed in this paper, ν e remained substantially smaller than ν t .Hence, rheological effects had no or only little impact on the simulation results.The bottom boundary condition to (2) can be formulated in terms of a friction velocity, u * , as: where τ b is bed shear stress, and ρ o = 1028 kg/m 3 is seawater density, noting that sediment concentrations of ~100 g/L increase the bulk density of the seawater-sediment mixture by ~6%, which is relatively small and can be ignored in (5).Using the law of the wall, the friction velocity is derived from the model velocity, u w + u b , at the first velocity grid point ( = ∆z/2) above the bed, yielding [38]: where κ = 0.4 is the von Kármán constant, and the roughness height K s parameterizes the thin region close to the bed not being resolved by the model.This study uses K s = 2 mm, which is close to the value (3 mm) used before [35].Note that ( 6), assuming the existence of a rough wall, is only valid if > z o where the roughness parameter z o is related to roughness height via z o = K s /30 [38].In our study, = 0.5 mm clearly exceeds z o = 2/30 mm ≈ 0.067 mm.The use of ( 6) also implies that the simulation results depend on the grid spacing of the bottom-nearest grid cell.Hence, when using a different grid spacing the erosion-rate constant e o in (12) needs to be re-adjusted in order to reproduce the laboratory mud suspensions.In the end, calibrated predictions using a different vertical resolution are very similar to those presented here (results not shown).

1DV Model: Suspended Sediment Equation
For simplicity, only a single sediment fraction of silt is considered in this study.The Reynolds-averaged equations for suspended sediment mass for mud-laden flow on a plane seafloor can then be written as: where C is suspended-sediment concentration, ν s is vertical sediment diffusivity, and w s is the settling velocity of floccules.Sediment diffusivity is derived from: where the turbulent Schmidt number is set to an often used value of σ s = 0.7.However, reported values of σ s are within a wide range from 0.25 to 2 [39], which the present study addresses in a sequence of sensitivity experiments.The settling flux, w s C, in ( 7) is subject to hindered settling.This behavior can be parameterized as e.g., [1]: where w so is the maximum settling velocity of mud floccules and C* is the value of C for which the settling flux is zero.The value of C* can be derived from the observed maximum of the settling-flux curve (9) that corresponds to C m = C*/6.Note that C* is different from the gelling concentration, C gel , i.e., the concentration at which floccules form a space-filling network.Settling experiments with cohesive sediment suggest that C gel is substantially smaller than C* with values of C m /C gel = 0.4 ± 0.1 [40].This relation has an upper bound of C gel ≈ C*/2, which is taken as gelling point in the following.
Permeability effects affiliated with the onset of consolidation come into play at higher sediment concentrations.These effects can be incorporated in the settling-flux parameterization via the refinement [40]: for C > C * /χ (10) where χ > 1 is an empirical parameter.We use χ = 3, which implies that permeability effects come into play before the gelling point at C*/2 is reached.Previous high-resolution studies e.g., [35] ignored permeability effects, which in the absence of other processes imply a small residual settling flux in high-density suspensions near the gelling point.

1DV Model: Sediment Flux across the Seafloor
The sediment flux across the seafloor can be formulated as: where D and E are rates of deposition and erosion, respectively.The erosion flux E can be calculated from the classical Ariathurai-Partheniades formula e.g., [41]: where e o is the erosion-rate constant, τ b is bed shear stress, and τ s is bed shear strength.Previous flume experiments have shown that e o is highly variable as it depends on the bulk density (or water content) of the sediments, particle size distribution as well as mean particle size, mineralogy, organic content, and amounts and sizes of gas bubbles see [42].In this study, the value of e o was empirically derived from calibration against the laboratory findings [22,23] with the target to capture the lower end of the observed suspended sediment concentrations.Erosion only occurs for τ b > τ s .Based on early laboratory findings [43], the bed shear strength τ s varies with erosion depth d (defined in relation ( 14)) according to: where τ so is a trigger value, and the shear-strength parameter d o refers to the depth of erosion where τ s = 2τ so .We assume that wave loading has already reworked and liquefied the seabed, so that sediment erosion is supported by relatively small values of τ so of 0.05-0.1 Pa, as suggested from laboratory experiments [2] and field data [21].On the other hand, the shear-strength parameter d o is chosen from a range of 0.1 to 0.25 mm, which characterizes a relatively "fresh" mud-dominated seabed of consolidation times <10 days [44].Additional interactive wave loading effects that could possible lower τ s during the process are ignored.Biological processes can also strongly modulate bed shear strength [45], and these are also ignored here.Equations ( 12) and ( 13) represent so-called Type 1 erosion [45] for cohesive sediment in which the bed shear strength is parameterized as a function of the amount of eroded sediment mass.Accordingly, erosion depth d (per unit surface area) is related to eroded sediment mass M according to: where ρ s,bed = 1150 kg/m 3 approximates the density of "soft" consolidated bed sediment corresponding to a sediment concentration of ~200 kg/m 3 .
Deposition of cohesive sediment is a continuous process independent of bed shear stress [46].Hence, the deposition flux, D, is of the same form as the settling flux (10); that is, where C b is suspended sediment concentration in the bottom-nearest grid cell of the model.

1DV Model: Turbulence Closure
Vertical eddy viscosity, ν t , is derived from a k-ε turbulence model that accounts for effects due to suspended sediment [34].Vertical eddy viscosity is calculated from [47]: where l * is a turbulent mixing length, k is the kinetic energy of turbulent velocity fluctuations, and ν min = 10 −10 m 2 s −1 is a very small dummy value ensuring that the model's eddy kinetic energy production never fully vanishes in the presence of vertical shear flow.The turbulent mixing length in ( 16) is given by: where c µ = 0.09 is an empirical constant, ε is dissipation of turbulent kinetic energy, and l max is an upper bound of l * , being defined as distance from the seafloor (which is a bound for the size of the largest vortices).Turbulent kinetic energy, k, and dissipation rate, ε, are estimated by their balance equations [35]: where the diffusion parameters are given by ν k = ν m + ν t /σ k and ν ε = ν m + ν t /σ ε , P k denotes the production term, P b and S b represent turbulence damping via density stratification of ambient seawater (ignored here) and suspended sediment, respectively, and S ε describes turbulence damping via viscous drag effects of individual sediment particles.Turbulence production is given by: where u = u w + u .Sediment-induced stratification effects are related to: Dissipation in carrier fluid turbulence caused by viscous drag of particles in the fluctuating field can be formulated as [35]: where T p is the viscous particle response time, and T L is the turbulent eddy timescale.For relatively small mud floccules, studied in this work, the particle response time is generally much shorter than the turbulent eddy timescale; that is, T p << T L = k/(6ε) [35].In this situation, ( 22) can be approximated by: which constitutes a concentration-dependent dissipation term.The magnitude of this term starts to exceed that of the conventional dissipation term ε in (19) when C > 90 g/L.
The numerical parameters in (18) and ( 19) are set to their standard values (σ k = 1.0, σ ε = 1.3, c ε1 = 1.44, and c ε2 = 1.92) except for c ε3 , which is the weighting factor of the dissipation terms in (19).Consistent with previous studies that suggested a large value range of this parameter [48], here we vary the value of c ε3 between 0.6 and 3 in a sequence of sensitivity experiments.
Zero-gradient conditions are used for k at both the surface and bottom of the water column.A zero-gradient condition is applied for ε at the surface of the water column.The bottom boundary condition for ε is given by [49]: with the same parameter settings as specified for ( 5) and (6).Note that condition ( 24) is grid-spacing dependent, but the predicted turbulence characteristics do not dramatically change when the vertical grid spacing is varied within a reasonable range (results not shown).

Experimental Design
This study considers two different approaches.This first approach assumes that the fluid mud layer maintains a constant thickness of h = 4 cm (which is consistent with laboratory findings, see [23]) and a vertically uniform suspended-sediment concentration, C = C b .The further condition of stationarity (i.e., the fluid mud layer does not move to other locations of different seabed properties) then implies M = d ρ s = h C b , which allows for the formulation of a single equation that returns instant sediment fluxes, E − D, across the seabed as a function of bed shear stress, τ b , and concentration, C = C b , of the fluid mud layer.The resulting equation, which follows from Equations ( 12)-( 15), reads: where the bed shear strength is given by: Note that our previous study [12] analyzed the equilibrium values of C b , i.e., the roots of (25), in the absence of permeability effects.
In the second approach, the sediment flux across the seafloor derived from ( 12) and ( 15) is interactively coupled with the 1DV model and properties of the mobile fluid mud layer are predicted by this model.We use a standard model configuration (Table 1), which has been calibrated against laboratory findings [22,23].The model is hereby run for orbital wave speeds (u a in Equation ( 1)) in a range of 10 to 55 cm/s, and the only model parameter that is varied is C*, which modifies the sediment concentration at which the settling/deposition flux attains a maximum.While we have undertaken sensitivity studies for all model parameters, here we restrict the discussion to parameters that are expected to lead to marked changes in suspended sediment characteristics, such as variations of w so or C*, and parameters of high uncertainty, such as the turbulent Schmidt number, σ s , or the stratification parameter, c ε3 in the ε Equation (19).Parameters related to the bed shear stress (6) and erosion flux (12) are kept unchanged throughout this study.Note that the apparently high value of the erosion-rate constant of e o = 0.01 kg/m 2 /s in conjunction with the other parameter settings is essential to achieve reasonable agreement with the laboratory data.For academic purposes, we also present experiments in which either the turbulence damping effect by sediment stratification (21), or that by individual sediment particles (23), is set to zero.

Results
We present three major results.The first is an analysis of previous laboratory data [22,23] in more detail.The second is derived from the stand-alone analysis of sediment fluxes across the seafloor for a fluid mud layer of constant height and uniform sediment concentration.The third is an analysis of the predictions by the 1DV model.

Characteristics of Mud Suspensions in Laboratory Experiments
Variations of the maximum horizontal speed of the oscillatory flow in laboratory experiments [22,23] resulted in a large range of equilibrium near-bed suspended sediment concentrations, C b (Figure 3).From the results we can identify three different regimes, noting that mud suspensions exclusively formed when wave velocities exceeded a threshold of ~20 cm/s.In the first regime, the relationship between wave velocity and C b is approximately linear for wave velocities between 20 and 35 cm/s.The second regime consists of four experiments in which C b remains limited to a value of ~50 g/L, despite the increase in wave velocities from 40 to 55 cm/s.In the third regime, there are five data points covering a wide range of equilibrium concentrations of C b from 22 to 80 g/L despite only small variations of wave velocities from 37 to 41 cm/s.This large range of C b despite small variations of wave velocities is also seen in the field data [1], but not necessarily for the same reasons.It is plausible to hypothesize that this wide range of C b despite similar forcing conditions is related to either a modification of associated bed shear stresses via the creation of bed ripples, or variations in overall turbulence characteristics of the wave boundary layer.A closer inspection of the laboratory results, however, reveals that neither of these factors explains the observed variability (Figure 4).Instead, the large range of C b values concurs with a narrow range of bed friction velocities of ~7-8 cm/s (vertical line in Figure 4a) and similar Reynolds numbers of ~1 × 10 5 (vertical line in Figure 4b).Hence, the observed variations of C b are for similar wave velocities of ~40 cm/s are either due to unknown experimental biases, or a feature of variable turbulence levels within the fluid mud layer not captured by the definition of the Reynolds number for the wave boundary layer.On the other hand, there seems to be a second regime in which C b reaches a quasi-equilibrium value of ~50 g/L, despite a relatively large range of bed friction velocities from 4.5 to 10.5 cm/s (Figure 4a) and Reynolds numbers from 1.5 × 10 5 to 2.7 × 10 5 (Figure 4b).This feature is indicative of the situation whereby the mud suspension has reached the gelling point, so that it cannot take up more bed sediment.

Erosion-Deposition Characteristics of Suspended Mud
The following discussion is based on the use of (25) with parameter settings of the main configuration (Table 1), a thickness of the fluid mud layer of h = 4 cm, and either C* = 180 g/L, or C* = 400 g/L.Initially we assume there is no suspended sediment present.Consequently, bed erosion with τ b > τ s will trigger a net sediment flux across the bed (E − D) such that the suspended sediment concentration increases at a rate of: This will continue until the net sediment flux across the bed eventually vanishes.The final equilibrium situation corresponds to the first crossing of erosion-flux and deposition-flux curves (Figure 5).With C* = 180 g/L, the bed shear stress of τ b = 0.5 Pa leads to an equilibrium suspended sediment concentration of C b ≈ 15 g/L (Figure 5a).An increased bed shear stress of τ b = 0.8 Pa triggers an equilibrium at C b ≈ 45 g/L.Further increases in τ b imply that the erosion flux always exceeds the deposition flux irrespective of the level of suspended mud concentration.According to (26), the suspended mud concentration will then continue to increase at a timescale that strongly depends on the minimum value of (E − D) encountered during the process.This is the gelling ignition mechanism, see [12].To limit C b to the gelling concentration at C gel = C*/2, which is not accounted for in the model so far, the erosion rate needs to be constrained.Here we achieve this by multiplying the erosion flux (12) with the limiting function: (28) which assumes that erosion fluxes vanish in close proximity of the consolidating bed.
In the second case, C* = 400 g/L, the same sequence of forcing (τ b = 0.5, 0.8, 1.5 Pa) leads to equilibrium values of C b = 13, 24, 60 g/L, respectively (Figure 5b).The gelling ignition mechanism is not triggered in this case.Overall, mud concentrations are smaller compared to the first case with the same forcing, but gelling ignition would create a much higher concentration of ~200 g/L.While the required bed shear stresses may appear unusually high, oscillatory currents of short periods (such as wind waves) can induce very high instantaneous bed shear stresses in shallow water e.g., [8].

One-Dimensional Water-Column Model Simulations
Here we present the results of the 1DV model simulations.First, we present a sequence of model simulations using the main configuration of parameters (detailed in Table 1) in comparison with laboratory data [22,23].Then, we discuss a number of sensitivity studies with a focus on parameters that control the thickness of the fluid mud layer.

Variations of Wave Forcing and C*
The applied wave forcing leads to the establishment of equilibrium suspended sediment profiles in all experiments.Overall, the simulated near-bed sediment concentration values agree well with the range measured in the laboratory (Figure 6).In each case, the equilibrium value of C b corresponds either to the first root of the boundary equation for sediment flux across the bed or to the assumed gelling concentration.The onset of the gelling ignition mechanism can be identified by a sudden sharp increase of C b for a relatively small increase in u b (Figure 6).The equilibrium value of M ranges between 1.5 and 3 kg/m 2 (depends on C*).The lutocline height, expressed by h = M/C b , approaches a similar value of ~3.5 cm, irrespective of the value chosen for C* (Figure 7c).
Let us further analyze the equilibrium state for a wave forcing of u a = 40 cm/s and C* = 180 g/L.Despite establishment of an overall equilibrium, the suspended sediment concentration of the fluid mud layer oscillates marginally by ±0.35 g/L on the time scale of the wave forcing (Figure 8a).The near-bottom oscillatory flow closely follows the ambient wave forcing (Figure 8b).Due to frictional effects, however, it experiences an amplitude reduction by ~45% and a slight phase shift.Nevertheless, the variations of the ambient flow field are too rapid to cause significant frictional flow reductions near the seabed.As a consequence of this, the oscillating flow still creates a bed shear stress of <τ b > ≈ 1.16 Pa on average, oscillating between zero and a maximum value of ~2 Pa on a time scale of 3 s, i.e., half the wave period (Figure 8c).This average bed shear stress exceeds the bed shear strength τ s ≈ 0.44 Pa by far, but, on average, the deposition flux operates to balance the resultant erosion flux (Figure 8d).On the time scale of half the wave period, net erosion fluxes vary by ±10 g/(m 2 s), which trigger the oscillations in C b .

How Do Lutoclines Form Under Wave Forcing?
All experiments using the main parameter configuration lead to the creation of high-density mud suspension with a pronounced lutocline structure (Figure 9a), which is explored in more detail in the following.Previously, it has been suggested that lutoclines form due to either modified settling speeds [50] or diffusive fluxes [28,29].The first hypothesis relies on increased settling speeds at lower suspended sediment concentrations, which leads to a vertical convergence of the settling flux.The second hypothesis relies on turbulence damping due to sediment stratification effects.Here, a wave forcing of u a = 40 cm/s creates a lutocline at ~3.5 cm above the bed (Figure 9a).This is of the order of the thickness of the wave boundary layer (Figure 9b), which can be estimated from the location at which the vertical gradient of the RMS values of the streamwise velocity variations vanishes.Sediment concentrations are too low to induce noticeable rheological effects.Hence, turbulent effects play a dominant role in vertical sediment diffusion inside the fluid mud layer which is consistent with laboratory observations [22,23].To this end, the turbulent sediment diffusivity attains a maximum value of ν s ≈ ν t /σ s ≈ 7 × 10 −5 m 2 /s with a few mm from the seafloor and then gradually decreases to vanish at the lutocline (Figure 9c).
The first hypothesis suggesting a vertical sediment flux convergence at the lutocline [50] does not apply here, simply because there is no source of suspended sediment in the ambient fluid above the lutocline that could balance the settling-induced sediment loss.Hence, lutocline creation can only be explained via modification of vertical sediment diffusion, as further explored in the following.The vertical sediment diffusion term can be expressed by two components: which can be interpreted as (1) a conventional diffusion term associated with a constant diffusivity; and (2) a correction term associated with vertical variations in the diffusivity.Surprisingly, it turns out that the conventional diffusion term is negative in most of the fluid mud layer (Figure 10a).Hence it is the last term in (29), which has a pronounced maximum at the lutocline (Figure 10b), that operates to counteract the effect of sediment settling.This simulation result suggests that vertical sediment stratification is responsible for the formation of lutoclines, not through localized turbulence damping, but in combination with a vertical decrease of sediment diffusivity (see Figure 9c).With the parameter settings of the main configuration (Table 1), the simulated equilibrium thickness of the fluid mud layer, h, lies in a range of 3-5 cm (Figure 11).This order of magnitude of h agrees well with the laboratory findings, which is not surprising, given that this was a goal of the model calibration.The lutocline height increases linearly with u a according to the model results, but this relation is only partially supported by the laboratory data.[23], Table 1.
While variation of C* does not significantly influence h (Figure 9a), the following section outlines model parameters and associated processes that control the equilibrium thickness of wave-created fluid mud layers in the model simulations.Constant values of ν s and w s lead to classical exponential Rouse profiles of sediment suspensions of the form of [51]: where z* is distance from the bed, which implies an exponential length scale of ν s /w s .While the assumptions behind (30) are grossly oversimplified, it is reasonable to suggest that w so in (10) and the turbulent Schmidt number, σ s , in (8), are the prime suspects for influencing h.Surprisingly, when varied over one order of magnitude, neither of these two parameters induces significant changes in h with the model applied here (Figure 12a,b).For instance, a tenfold reduction of w so from 1 mm/s to 0.1 mm/s increases h from 3.5 to 4.7 m, which is by a factor of only 1.34.On the other hand, the value of the Schmidt number does not seem to have an influence on h at all.After exploring variations of all model parameters within reasonable ranges, we concluded that only a single parameter, namely c ε3 , which controls the relative magnitude of sediment stratification effects in the ε-equation ( 19), can influence h markedly.When selecting c ε3 from the previously reported range of values from 0.6 to 3 see [48] we yielded large variations of h from <1 cm to >8 cm (Figure 12c).
Finally, in order to explore the difference between fluid mud layers of pronounced lutoclines and exponentially shaped sediment profiles, we repeated the main experiment with u a = 40 cm/s and C* = 180 g/L for another two cases.Case 1 is a model run without particle-induced turbulence damping, i.e., S ε = 0 in ( 18) and ( 19).Case 2 ignores sediment stratification effects in the turbulence closure scheme, i.e., S b = 0, noting that the structure of the wave boundary layer remains unaffected by these modifications (Figure 13b).Case 1 only produces minor changes in the results.Case 2 leads to dramatic changes.In this case, turbulent eddy viscosity and, hence, sediment diffusivities are no longer constrained by stratification effects at the lutocline.Hence, turbulence can develop outside the wave boundary layer (Figure 13c) and enhanced vertical sediment diffusion leads to a gradual thickening of the suspended sediment layer (Figure 13a).Due to continued diffusive upward sediment fluxes, the near-bed suspended sediment concentration does not reach its equilibrium value, but continues to increase over the simulation (Figure 14a).So do the total eroded sediment mass, M, and the vertical length scale, h = M/C (Figure 14b,c).Hence, the damping of turbulence by sediment-induced density stratification via the term S b in ( 18) and ( 19) is fundamental in the creation of lutoclines, which is consistent with earlier suggestions [28,29].

Discussion
An important finding of this study is that wave-created mud suspensions attain one of two possible equilibrium suspended sediment concentrations, either the gelling concentration, C*/2, for high-enough bed shear stresses, or a value that lies below C*/6, corresponding to the maximum of the deposition-flux curve.The transition phase between these two values has been previously referred to as gelling ignition [12], which can be clearly identified in our advanced 1DV simulations.It is important to note that equilibrium states cannot always be interpreted as gelling concentrations, which makes the interpretations of measurements difficult.
Overall our model predictions agree with previous laboratory studies [22,23] and limited field data [1].The model parameter C*, which defines the maximum of the settling/deposition-flux curve, has a dominant influence on the equilibrium concentration, C b , of mud suspensions.Thus, large variations in C b for similar wave velocities can be explained by variations of C*.Previous laboratory studies have demonstrated that ambient turbulence can significantly shift C* and, hence, the highest possible suspended sediment concentration to higher values [18].In the present 1DV model, C* needs to be prescribed as an external parameter, and one of the key challenges for future studies is to link C* to turbulence properties within the fluid mud layer.While we cannot exclude experimental inconsistencies in the laboratory findings, the fundamental equations governing the 1DV model applied in this work suggest that variations in mud concentrations under similar wave forcing could be related to changes in C*, possibly due to changes of turbulence levels inside the mud layer.More laboratory experiments are required to verify this hypothesis.
Another key finding of our study is that the upward diffusive flux in wave-created fluid mud peaks at the lutocline.This result is in contrast to the "saturation hypothesis" [29], postulating that sediment-induced turbulence damping at the interface between the fluid mud layer and the overlying fluid creates a total collapse of the turbulent energy and, hence, the disappearance of vertical diffusive sediment fluxes.On the other hand, the findings presented here also indicate that the dissipation of turbulent energy via sediment-induced density stratification has a dominant effect on the vertical sediment diffusivity and, hence, the vertical structure of wave-created fluid mud suspensions.In the end it turns out that the interior of wave-created mud suspensions is highly turbulent while the lutocline prevents turbulence from radiating into the ambient fluid above the fluid mud layer (see Figure 13c, Case 1).This finding is consistent with previous interpretations of high-resolution turbulence measurements in wave-created mud suspensions [22].
While the validity of the k-ε turbulence closure may be questioned in the context of fluid mud layers, findings presented here confirm that its use within the 1DV modelling framework leads to reasonable predictions of lutocline structures.It should also be noted that high-resolution turbulence measurements [22] demonstrate that the sediment suspensions in the laboratory experiments discussed here were fully turbulent, which is consistent with the predicted eddy-viscosity profiles (see Figure 9c).
It should be pointed out that an additional treatment (28) was required in our model to prevent the unrealistic increase of the concentrations of mud suspensions beyond the gelling point.Hereby we assumed that erosion fluxes tend to vanish near the gelling point, which needs to be tested in future investigations.
One interesting feature that not included in our 1DV model is the diffusive leakage or detrainment of suspended sediment across the lutocline into the ambient fluid.This leakage is obvious in laboratory experiments [2,23] where suspended mud concentrations >1-2 g/L developed above the lutocline.Interfacial instabilities on the lutocline due to the Kelvin-Helmholtz instability mechanism e.g., [52], not accounted for in this study, may cause such a leakage-to be studied in more detail in the future.Ambient steady flows, not considered in the present study, are on their own not energetic enough to induce bed shear stresses >1 Pa required for the creation of fluid mud layers via the erosion mechanism.Hence, in the overall context of suspended sediment dynamics, ambient steady flows are more relevant as supply agents carrying fine sediment from rivers to the continental shelf.On the other hand, wave-created fluid mud suspensions may accelerate on a sloping seafloor to form a type of self-accelerating turbidity current [53,54] in the form of a self-moving mud slurry [19], which is worth further investigation.

Figure 2 .
Figure 2. Schematic of the gelling ignition mechanism, which is the transition from low to high suspended mud concentrations, after [12].Shown are typical distributions of erosion and deposition fluxes as a function of suspended sediment concentration C of a lutocline layer of constant thickness and uniform C.Under low bed shear stresses (solid erosion curve), net erosion (E-D > 0) occurs until C has reached the equilibrium concentration of C 1 , giving the lower limit of the "depositional hurdle", which prevents C from increasing any further.C 2 is the second (unstable) root of the erosion-deposition equation.C m refers to the maximum of the deposition flux.The depositional hurdle is overcome for increased bed stresses (dotted erosion curve), and C will increase unhindered until reaching the gelling point.

Figure 3 .
Figure 3. Equilibrium near-bed concentration (C b ) as a function of wave speed (u a ) measured in the laboratory [23], Table1.Also shown are field data (triangles)[1].Dashed lines indicate different regimes, as discussed in the text.

Figure 4 .
Figure 4. Equilibrium near-bed concentration (C b ) as a function of (a) bed friction velocity (u* gm ), derived from an empirical formula [8] taking into account for bed ripples; and (b) the Reynolds number for the wave boundary layer.All data are taken fromTable 1 of [23].Dashed lines indicate different regimes, as discussed in the text.

Figure 5 .
Figure 5. Erosion and deposition fluxes in (25) as a function of bed shear stress, τ b , plotted against C/C* for the parameter settings of the main configuration (see Table 1) and an assumed constant lutocline height of h = 4 cm, for (a) C* = 180 g/L and (b) C* = 400 g/L.Crossings between curves imply equilibrium states of zero fluxes.Arrows indicate the corresponding values of C/C*.Additional treatment is required in the regime highlighted by the circle in (a) to limit concentrations by the gelling concentration at C/C* ≈ 0.5.See text for more details.The equilibrium mud concentrations in (a) are 15, 45 and C*/2 = 90 g/L; and in (b) 13, 24 and 60 g/L.
Figure 5. Erosion and deposition fluxes in (25) as a function of bed shear stress, τ b , plotted against C/C* for the parameter settings of the main configuration (see Table 1) and an assumed constant lutocline height of h = 4 cm, for (a) C* = 180 g/L and (b) C* = 400 g/L.Crossings between curves imply equilibrium states of zero fluxes.Arrows indicate the corresponding values of C/C*.Additional treatment is required in the regime highlighted by the circle in (a) to limit concentrations by the gelling concentration at C/C* ≈ 0.5.See text for more details.The equilibrium mud concentrations in (a) are 15, 45 and C*/2 = 90 g/L; and in (b) 13, 24 and 60 g/L.
The associated value of C b coincides with the maximum of the deposition-flux curve, i.e., C = C*/6, which depends on the setting of C*.To this end, we can match the two instances of the highest concentrations >70 g/L observed in the laboratory with the choice of C* = 180 g/L.Laboratory cases in which mud concentrations stagnate at a level ~50 g/L can be largely reproduced with the choice of C* = 90 g/L.On the other hand, the setting of C* = 400 g/L gives a curve that agrees well as a lower bound of the laboratory results.In general, variation of C* explains the observed large spread of C b values seen for a wave forcing of ~40 cm/s.Note that our model reproduces high fluid mud concentrations >120 g/L, observed in the field [1], with the setting of C* = 240 g/L.

Figure 6 .
Figure 6.Same as Figure 3, but including model-simulated equilibrium near-bottom suspended sediment concentrations, C b , from 4 sequences of numerical experiments.

Figure 7 .
Figure 7. Main parameter configuration with a wave forcing of u b = 40 cm/s.Evolution of (a) near-bed suspended-sediment concentration, C b ; (b) total eroded sediment mass, M; and (c) lutocline height, h = M/C for C* = 180 g/L (dashed black line) and C* = 400 g/L (solid red line).

Figure 8 .
Figure 8. Main parameter configuration with u a = 40 cm/s and C* = 180 g/L.Shown are sequences of 12 s in duration starting after 10 min of simulation of (a) near-bed suspended-sediment concentration; C b ; (b) near-bed horizontal velocity, u w + u' b ; (c) bed shear stress, τ b , and bed shear strength, τ s ; and (d) net sediment flux across the seabed, E − D. The average bed shear stress is <τ b > ≈ 1.16 Pa.The average bed shear strength is <τ s > ≈ 0.44 Pa.

Figure 9 .
Figure 9. Main configuration with u a = 40 cm/s and C* = 180 g/L (dashed lines) or C* = 400 g/L (solid lines).Equilibrium vertical profiles of (a) suspended sediment concentration, C; (b) horizontal velocity, u w + u', averaged over a wave period, and (c) vertical eddy viscosity, ν t , noting that effective sediment viscosity remains insignificantly small, ν e < 0.1 × 10 −4 m 2 /s throughout the fluid mud layer.z* is distance from the seabed.This figure shows the dynamical properties of the lutocline layer.

Figure 10 .
Figure 10.Main configuration with u a = 40 cm/s and C* = 180 g/L.Shown are the equilibrium profiles of the two terms on the right-hand side of (29), making up the vertical sediment diffusion term in(7).The insert in (b) shows details for the lower 3 cm of the water column for comparison with (a).

Figure 11 .
Figure 11.Lutocline height, based on the location of the C = 10 g/L level, as a function of wave forcing, u a , for the main parameter configuration (diamonds, setting of C* is irrelevant) and laboratory experiments (circles)[23], Table1.

Figure 12 .
Figure 12.Sensitivity studies.Lutocline height, based on h = M/C b , for a wave forcing of u a = 40 cm/s and C* = 180 h/L as a function of (a) maximum settling speed, w so , in (10) and (15); (b) the Schmidt number, σ s , which appears in (8) and (21); and (c) the stratification parameter c ε3 in(19).In all cases, C b attains an equilibrium value of 84 ± 4 g/L.The value of c ε3 has a dominant effect on lutocline height.

Figure 13 .
Figure 13.Sensitivity study.Vertical profiles of (a) suspended sediment concentration, C; (b) profile of RMS values of horizontal velocity variations, u w + u , and (c) vertical eddy viscosity, ν t , after 30 min of simulation.z* is distance from the seabed.Shown are case 1 (solid black lines; S ε = 0) and case 2 (dashed red lines; S b = 0).Sediment stratification effects are essential for lutocline development.

Figure 14 .
Figure 14.Same as Figure 7 with C* = 180 g/L, but showing results for case 1 (solid black lines; S ε = 0) and case 2 (dashed red lines; S b = 0).Sediment stratification effects have a dominant impact on lutocline height.

Table 1 .
Key parameters and values.