Critical Shear Stress for Erosion of Sand-Mud Mixtures and Pure Mud

The erosion threshold of sand-mud mixtures is investigated by analyzing the momentum balance of a sand particle or a mud parcel in the mixture bed surface, and a formula for the critical shear stress of sand-mud mixtures is developed, which also applies for pure sand and mud. The developed formula suggests that the variation of the critical shear stress of sand-mud mixtures over mud content is mainly caused by the varying dry bulk density of the mud component in the mixture. The developed formula reproduces well the variation of the critical shear stress of sand-mud mixtures over mud content and can predict the critical shear stress of both sand-mud mixtures and pure mud in the process of consolidation. The developed formula promises to be convenient for application by relating the critical shear stress to mud content and the dry bulk density of sediment.


INTRODUCTION
The erosion of sediment is one of the controlling processes of sediment dynamics in aquatic systems. It concerns various domains, including geomorphology, pollutant transport, dredging activities, scour around structures, landward retreat of shorelines, bank failure, etc. (Sanford and Maa, 2001;Mostafa et al., 2008;Chen et al., 2017;Kurdistani et al., 2019;Zheng et al., 2020;Li et al., 2021). Assessment of the erodibility of sediment has therefore been of interest to numerous scientists and engineers. The particle size is one of the most critical factors affecting the erodibility of sediment. According to the particle size, sediments can be classified into gravel, sand, silt, and clay. Coarse-grained sediments, such as coarse silt, sand and gravel, are usually non-cohesive. Sediments comprised primarily of fine silt and clay exhibit cohesion, and therefore are often referred to as cohesive sediments or mud (Mehta and Partheniades, 1982;Van Rijn, 1993). The physicochemical properties of non-cohesive sediment and cohesive sediment are strikingly different. As a result, the erodibilities of non-cohesive and cohesive sediments have been separately studied in the past several decades (Van Ledden et al., 2004;Le Hir et al., 2011).
Non-cohesive sediment is eroded as individual particles. The erosion resistance of non-cohesive sediment is mainly provided by the submerged weight of the sediment, which leads to express the erosion threshold as a function of the size, shape and density of the particles (Miller et al., 1977;Chiew and Parker, 1994). For cohesive sediment, the interparticle attractive force coming from electrochemical effects is much more significant than the gravitational force and plays a dominant role in the erosional behavior and property (Righetti and Lucarelli, 2007). Experimental studies have shown that cohesive sediment is eroded as aggregates (i.e., surface erosion) at low bed shear stresses and as lumps or chunks of bed material (i.e., mass/bulk erosion) at high bed shear stresses (Partheniades, 1965;Mehta, 1989;Sanford and Maa, 2001;Winterwerp et al., 2012). The erosion threshold of cohesive sediment has been found to increase with increasing degree of compactness of sediment when other physicochemical properties of the sediment are controlled. Some empirical relationships were also developed between the critical shear stress and physical properties of cohesive sediments, including wet or dry bulk density (Tang, 1963;Owen, 1970;Dou, 2000;Amos et al., 2004;Xu et al., 2015;Zuo et al., 2017), water content (Jacobs et al., 2011), solid volume fraction (Kusuda et al., 1984;Mehta and Parchure, 2000), solid/void volume ratio (Wu et al., 2017) and yield stress . As the erodibility of cohesive sediment is also affected by factors like mineralogy, pore water chemistry, temperature, organic content, biological actions, etc., those relationships between the erosion threshold and the physical properties are often modified by these influences. Summaries on the theories and models for predicting the erosion threshold and rate of non-cohesive and cohesive sediments could be found in the studies of Zhu et al. (2008); Le ), and Chen et al. (2018. In many natural environments, e.g., estuary, delta, mangrove forests, muddy and silty coasts, non-cohesive sediment and cohesive sediment are not completely separated and often occur as sand-mud mixtures (Carniello et al., 2012;Mehta and Letter, 2013). Compared with pure sand and pure mud, the erodibility of sand-mud mixtures has not been well-understood. A few experiments have been conducted to investigate the erodibility of sand-mud mixtures in the past decades (Torfs, 1995;Panagiotopoulos et al., 1997;Sharif, 2003;Le Hir et al., 2008;Jacobs et al., 2011;Ye et al., 2011;Smith et al., 2015;Van Rijn, 2020). These experiments have shown that sand-mud mixtures behave like cohesive sediments if they contain enough cohesive particles, otherwise like non-cohesive sediments (Torfs, 1995;Panagiotopoulos et al., 1997). The critical shear stress of sand-mud mixtures has been found significantly correlated to mud content. For non-cohesive mixtures, the critical shear stress generally increases with increasing mud content (Torfs, 1995;Le Hir et al., 2008;Smith et al., 2015). In some cases, a slight decrease in the sand threshold was sometimes observed for very low mud fractions (Torfs et al., 2000;Barry et al., 2006). Barry et al. (2006) attributed this decrease to viscous lubrication induced by clay-sized particles in the bed pore fluid. For cohesive mixtures, the critical shear stress increases monotonously with increasing mud content or firstly increases with mud content up to an optimum mud content then decreases slowly with further increasing mud content until 100% mud (Nalluri and Alvarez, 1992;Sharif, 2003;Ye et al., 2011;Smith et al., 2015;Van Rijn, 2020).
Those experiments of sand-mud mixtures have shown that the mixtures behave differently from pure sand and pure mud. Applying the formulae developed for pure sand and pure mud to sand-mud mixtures would induce massive errors. Van Ledden (2003) developed an empirical formula for the critical shear stress of sand-mud mixtures, which is expressed as a function of the mud content, the critical mud content, the critical shear stresses of pure sand and pure mud: where τ cr is the critical shear stress of sand-mud mixtures; p m is mud content; p mcr is the critical mud content (about 10-15%), below which the mixture is cohesionless and above which the mixture exhibits cohesion; τ crs and τ crm are the critical shear stresses of pure sand and pure mud, respectively; β v is an empirical coefficient between 0.75 and 1.25. Ahmad et al. (2011) developed a similar formula for the critical shear stress of sand-mud mixtures. Compared with the formula of Van Ledden (2003), their formula is simpler: where p s is sand content (i.e., the fraction of sand in the mixture); and ζ is an empirical coefficient between 0.1 and 0.2. The formulae of Van Ledden (2003) and Ahmad et al. (2011) involve the critical shear stress of pure mud. However, both Van Ledden (2003) and Ahmad et al. (2011) did not give a solution to determine the critical shear stress of pure mud for a given sandmud mixture, which limits the application of their formulae in practice. Recently, Wu et al. (2017) developed a formula for the critical shear stress of sand-mud mixtures by analyzing the force balance on a sediment parcel beginning to erode from the bed surface: where τ crL is the critical shear stress for mixtures with low mud contents, τ crL = τ crs + 1.25 (τ crmc − τ crs ) min p m , 0.05 ; τ crmc is the critical shear stress for erosion of mud corresponding to the porosity of the mud component in a sediment mixture; α and β w are two coefficients. The value of α is related to the sand median diameter d s : α = 0.42 exp −3380d s . The value of β w is set as constant 1.2. Wu et al. (2017) assumed τ crmc could be calculated by the formula for the critical shear stress of pure mud. Based on their collected data of pure mud, they developed an empirical formula for the critical shear stress of pure mud: where r is the solid/void volume ratio, r = (1 − φ m ) /φ m , with φ m being the mud porosity. As the coefficient 10.29 in Eq. (4) may vary for different muds, Wu et al. (2017) suggested τ crmc could be calculated by: where τ crc is a known critical shear stress for reference, which can be the critical shear stress in the case of 100% mud or another high mud content if the pure mud is not tested; r c is the solid/void volume ratio of the mixture/mud corresponding to the reference critical shear stress. Chen et al. (2018) also developed a formula for the critical shear stress of sand-mud mixtures based on the analysis of the balance of forces on a single particle of the mixture: where d and ρ s are the representative size and density of primary particles of the mixture, respectively; θ cr0 is the critical Shields parameter of non-cohesive sediment which is of a diameter d and particle density ρ s ; ρ is density of water; d m and ρ pm are the diameter and density of primary particles of the mud component of the mixture, respectively; ρ d is the dry bulk density of the mixture; ρ ps is the density of sand particles; ρ sdm is the stable dry bulk density of the mud component (which is the dry bulk density of the mud component when it gets fully consolidated); θ cr1 and m are coefficients, θ cr1 = 6.20 × 10 −8 m 3 s −2 and m = 1.55. The formula of Chen et al. (2018) applies not only to sand-mud mixtures but also to pure sand and mud. However, the stable dry bulk density of the mud component involved in the formulation is difficult to determine for a given sand-mud mixture, limiting the application of their formula.
The present study aims to develop a formula for the critical shear stress of sand-mud mixtures. The developed formula promises to be simple and easily applied in practice. It also should cover the full range of the mud content, i.e., it applies for pure sand, pure mud and sand-mud mixtures. The formula development, testing and discussion are described in the following sections.

Particle Cohesion
The cohesion between the fine-grained particles arises from electrochemical effects and is often modulated by biochemical factors (e.g., mucopolysaccharide binding) (Mehta and Lee, 1994). The biochemical effect is not considered in this study, leaving a focus on the particle cohesion from electrochemical actions. The fine-grained particles usually carry net negative charges which attract cations in water to form a double-layer water film coating the particles. The coating water film cannot transfer hydrostatic pressure. Therefore, the water pressure would induce an additional force acting on the overlapping area when the water films coating the neighboring particles overlap. This additional force induced by water pressure has been verified by the experiment with cross-quartz fibers (Deriagin and Malkin, 1950). Some researchers have considered the additional force induced by water pressure as one of the origins of the cohesive force between fine-grained particles (e.g., Dou, 1962Dou, , 2000Han, 1982;Zhang, 2012). Since there is a fundamental difference between the force induced by water pressure and the cohesive force arising from electrochemical actions, we prefer to consider them two different forces. The additional force induced by water pressure is not taken into account in this study because most of the existing erosion tests of sand-mud mixtures and pure mud were conducted in small-depth water flumes.
The van der Waals attraction has been generally believed to be responsible for particle cohesion from electrochemical effects in cohesive sediments (Han, 1982;Lick et al., 2004;Righetti and Lucarelli, 2007;Ternat et al., 2008). For two spherical particles of equal diameter, the van der Waals force f c between them is given by (Hamaker, 1937): where d m is the diameter of the cohesive particles; l is the separation distance between two particles (i.e., the smallest distance between the surfaces of the particles); A h is the Hamaker constant which reflects the strength of the van der Waals force, with its value typically between 10 −19 and 10 −21 J. The Hamaker constant A h is a function of the interacting particles and the intervening medium (Mehta, 2014). Therefore, its value is site-specific as the mineral compositions and pore water environments of cohesive sediments from different sites are usually different. The van der Waals force is a short-range force with its effective acting range typically within ∼0.1 µm (Hoath, 2016). As the fine cohesive particles usually form loosely structures called aggregates or flocs, the separation distance between neighboring particles could be far larger than the effective acting range of the van der Waals force, i.e., the van der Waals force is not always effective between neighboring particles. Studies have shown that the van der Waals force is significant between two contacted cohesive particles and negligible between two particles not contacted directly (Han, 1982). Figure 1 shows the distinction between a separation distance between two contacted particles (l ) and a separation distance between two neighboring particles that are not contacted with each other (i.e., s − d m in which s is the distance between two neighboring particles). The separation distance between two contacted particles is generally on the same order of magnitude as the thickness of the water films coating the fine particles. The average separation distance s − d m between neighboring particles is usually on the same order of magnitude as the particle diameter. During the consolidation of cohesive sediment, the particle packing becomes dense and the cohesive force gets enhanced. Therefore, the average separation distance between two contacted particles is positively correlated to the average separation distance between two neighboring particles. And the two average separation distances decrease with increasing the compactness degree of the sediment. As a first approximation, it is assumed that the dimensionless average separation distance l /δ (where δ is the thickness of the water film coating the cohesive particles) between two contacted particles is proportional to the dimensionless average separation distance s − d m /d m between neighboring particles, i.e., where η is a coefficient. The average (center-to-center) distance between neighboring particles s could be estimated from the dry bulk density of the mud component (Ternat et al., 2008): Substituting Eqs. (8,9) into Eq. (7), the van der Waals force between two contacted cohesive particles could be estimated by:

Incipient Motion Analysis of Sand-Mud Mixtures
The erosion behavior of sand-mud mixtures has been found related to the network structure of mixtures (Van Ledden et al., 2004;Jacobs et al., 2011;Wu et al., 2017). There are two typical network structures of sand-mud mixtures. When the mud content is low, the sand particles contact each other and form a skeleton, with the mud particles filling the voids formed between the sand grains (Figure 2A). This kind of mixture (referring to as sand-dominated mixtures hereafter) behaves as non-cohesive or less cohesive (Torfs, 1995;Van Ledden et al., 2004;Wu et al., 2017). The erosion process of these mixtures is governed by the sand component and occurs in the form of the detachment of particles of sand and flocs of mud. When the mud content is high, there are sufficient mud particles to prevent grain-tograin contact of the sand particles, and consequently, the sand particles lose contact with each other and "float" in the mud matrix ( Figure 2B). This kind of mixture (called mud-dominated mixtures hereafter) behaves as cohesive and the typical two erosion modes of cohesive sediment can occur (Torfs, 1995;Smith et al., 2015;Wu et al., 2017).
Consider a flat horizontal sand-dominated mixture bed and a flat horizontal mud-dominated mixture bed exposed to unidirectional flow. The erosion processes of the sanddominated mixtures and the mud-dominated mixtures are, respectively, governed by the sand component and the mud component (Torfs, 1995;Van Ledden et al., 2004;Wu et al., 2017). Here we take the detachment of a sand particle in the bed surface as the incipient motion criterion of the sanddominated mixtures and the entrainment of a cohesive sediment parcel (which may be a cohesive sediment floc or aggregate) in the bed surface as the incipient motion criterion of the muddominated mixtures.
The sand particles in sand-dominated mixtures have been found usually coated by a thin layer of the mud particles (Revil et al., 2002;Duteil et al., 2020;Van Rijn, 2020;Worden et al., 2020). The phenomenon is often documented as grain coating and ascribed to the adhesive force (the attraction/bonding between two particles of different media), diagenesis, biological actions, etc. (Duteil et al., 2020;Worden et al., 2020). Here, we assume that the sand particles in the sand-dominated mixtures are coated by a thin layer of the mud particles, and when a sand particle is disrupted from the bed surface, the coating layer is entrained together with the sand particle (Figure 3). This assumption has an advantage that in the analysis of the incipient motion of the sand particle, the adhesive forces between the sand particle and mud particles could be left out of consideration, which makes it easier to study the incipient motions of the sanddominated mixtures and the mud-dominated mixtures in the same theoretical framework. Figure 3 shows the incipient motion condition of a sand particle (or a cohesive sediment parcel) belonging to a sanddominated mixture (or a mud-dominated mixture). The forces acting on the sand particle or the mud parcel include the drag force F d , the lift force F l , the effective gravitational force G and the additional force F c which is the resultant force of the interparticle attractive forces (the van der Waals forces) between FIGURE 3 | Forces acting on a sand particle (or a mud parcel) in the surface of a sand-dominated mixture (a mud-dominated mixture). For the sand-dominated mixture, regard the blue circle as a sand particle and the red circle as the coating layer; for the mud-dominated mixture, regard the red circle as a mud parcel. the mud particles in the surface of the parcel or the coating layer of the sand particle and those surrounding mud particles. Considering the entrainments of the sand particle and the mud parcel are usually completed within a very short period, both the sand particle and mud parcel are assumed rigid bodies. The momentum balance equation for the critical condition of the incipient motions of the sand particle and the mud parcel is given by: where d r is a representative diameter, which refers to the diameter of sand particles for sand-dominated mixtures or the average diameter of mud parcels for mud-dominated mixtures, i.e., d r = d s for p m ≤ p mcr and d r = d a for p m > p mcr , with d s being the diameter of the sand particle and d a being the diameter of the mud parcel; k 1 d r , k 2 d r , k 3 d r , and k 4 d r are the moment arms of the drag force F d , lift force F l , submerged weight G and resultant cohesive force F c , respectively, with k 1 , k 2 , k 3 , and k 4 set as the proportionality coefficients. The drag and lift forces are given by (Dou, 2000;Righetti and Lucarelli, 2007;Vollmer and Kleinhans, 2007;Wu et al., 2017), where C d and C l are drag and lift coefficients, respectively; ρ is the density of water; u b is the flow velocity near the bed surface and α 1 is the area shape factor of the sand particle or the mud parcel.
The effective gravitational force G is given by G = α 2 ρ pr − ρ gd 3 r , where α 2 is a volumetric shape coefficient; g is the gravitational acceleration and ρ pr is the density of the sand particle or the mud parcel, i.e., ρ pr = ρ ps for p m ≤ p mcr and ρ pr = ρ pa for p m > p mcr with ρ pa being the density of the mud parcel. For the sand-dominated mixtures, the effective weights of the mud particles in the thin coating layer are negligible due to the thin layer and the small size of the mud particles.
The resultant force F c could be obtained by integrating those cohesive forces coming from the mud particles surrounding the motion-initiating sand particle or mud parcel: F c = k 5 nc n f c , where f c is the cohesive force (the van der Waals force) between two contacted mud particles; n is the number of mud particles coating the motion-initiating sand particle or the number of mud particles in the buried surface of the motion-initiating mud parcel; c n is the coordination number, i.e., the average number of the contacted particles of a mud particle; k 5 is a coefficient. According to Meissner et al. (1964), c n is a function of the volume fraction of solids: where ρ dm ρ pm denotes the volume fraction of solids of the mud component. n could be calculated by: r denotes the buried surface area of the sand particle or the mud parcel in which η is the relative protruding height of the sand particle or mud parcel into the flow from the bed surface; and N is the number of mud particles per unit area of the buried surface of the sand particle or the mud parcel. Assuming the sizes of the sand particle and the mud parcel are far larger than the size of the mud particles, N could be estimated by: where d m and ρ pm are the diameter and density of the mud particles, respectively; s is the average distance between neighboring particles and ρ dm is the dry bulk density of the mud component of the mixture. Considering Eqs. (9, 13), n is given by: Considering Eqs. (10,12,14), the resultant force F c is given by: Substituting the expressions for F d , F l , G, and F c into Eq. (11), the near-bed flow velocity for the incipient motion of sediment, u b,cr , is given by: Assuming the log law of velocity is valid near the bed, the near-bed flow velocity could be calculated by (Wu et al., 2017;Zhang and Zou, 2019;Li et al., 2020): where z b is the elevation where the drag force acts on the sediment particle or parcel; χ is the correction factor of Einstein (Einstein, 1950) and k s is the equivalent roughness height. Both z b and k s are related to the size of the sediment particle or parcel. Substituting Eq. (17) into Eq. (16) and considering u * = √ τ b /ρ in which τ b is the bed shear stress, the critical shear stress of sandmud mixtures τ cr could be obtained: According to Eq. (18), the critical shear stress of sand-mud mixtures consists of two parts which are, respectively, contributed by the effective weight of the sand particle or the mud parcel (corresponding to the first term on the right side of the equation) and the cohesive strength from the mud component (corresponding to the second term). By ignoring the cohesion between particles [i.e., ignoring the second term in the right side of Eq. (18)], Eq. (18) is reduced to: The left hand side of Eq. (19) is the critical Shields parameter. This yields where θ cr0 d r * is the critical Shields parameter of non-cohesive sediment of a particle diameter d r and d r * is the dimensionless diameter, defined as d r * = d r ρ pr ρ−1 g/υ 2 1/3 . θ cr0 d r * could be calculated from the Shields diagram or by the formula of Soulsby and Whitehouse (1997): Substituting Eq. (20) into Eq. (18), the critical shear stress of sand-mud mixtures is given by: The erosion threshold of the mud-dominated mixtures mainly depends from the cohesive strength of the mud component. Therefore, the first term on the right side of Eq. (22) could be negligible for mud-dominated mixtures. As a result, Eq. (22) can be further written as where p m is mud content and p mcr is the critical mud content beyond which the mixture is a mud-dominated mixture or pure mud and below which the mixture is a sand-dominated mixture or pure sand. Equation (23) is the formula we develop for predicting the critical shear stress of sand-mud mixtures. It covers the full range of mud content from 0 to 100%, i.e., it also applies for pure sand and pure mud. When the mud content is 0%, Eq. (23) is reduced into the widely used formula for non-cohesive sediment: τ cr = θ cr0 d s * ρ ps − ρ gd s . When the mud content is higher than 0 but lower than the critical mud content, i.e., 0 < p m ≤ p mcr , the mixture occurs in a sand-dominated structure and its critical shear stress is a function of the diameter of the sand particles, the diameter of the mud particles and the dry bulk density of the mud component. When the mud content is beyond the critical mud content, i.e., p m > p mcr , the mixture occurs in a muddominated structure or as pure mud. For this case, the critical shear stress is a function of the diameter of the mud particles and the dry bulk density of the mud component (or of the pure mud).
When applying Eq. (23) to a specific mixture, the dry bulk density of the mud component in the mixture needs to be determined first. As the mud particles fill in the space between the sand particles, the dry bulk density of the mud component in a mixture could be estimated by ρ dm = ρ d p m (1 − ϕ s ) , where ϕ s is the volume fraction of sand particles. Considering ϕ s = ρ d p s ρ ps , the dry bulk density of the mud component in a mixture is then calculated by: Equation (24) shows the dry bulk density of the mud component is 0 for pure sand and equal to the dry bulk density of the mud when the mud content is 100%. There are two coefficients in Eq. (23): p mcr and A. According to the existing experiments of sand-mud mixtures, the erosion mode of mixtures changes from non-cohesive to cohesive at a mud content of approximately 5-15%. As some sand-dominated mixtures could also exhibit weak cohesive behavior, the upper value, 15%, is used for the critical mud content in this study, i.e., p mcr = 15% . The value of A not only reflects the cohesion strength of the mud component, but also is related to the roughness of the mixture bed surface. Currently, it is difficult to determine the value of A as some coefficients in the expression for A are usually unknown, e.g., A h , η , k 4 , and k 5 . Therefore, A is treated as an empirical coefficient that will be determined by the measured data of erosion thresholds.

VALIDITY OF THE DEVELOPED FORMULA IN SAND-MUD MIXTURES
Data of eight groups of sand-mud mixtures collected from four experiments conducted by Torfs (1995); Sharif (2003), Jacobs et al. (2011), andSmith et al. (2015) are used to test the validity of Eq. (23) in sand-mud mixtures. The synopses of the collected data are listed in Table 1. Eq. (23) is applied to each group of data in this section. A best-fit value of A is adopted for each dataset, which was obtained by the regression analysis of the data. We also test the formulae of Van Ledden (2003) Experimental Data of Torfs (1995) Torfs (1995) conducted a series of erosion experiments of artificial mixtures of sand and clay mineral in a straight flume with a rectangular cross-section. The mixtures were prepared by mixing clay mineral with fine sand and were kept at a constant bulk density around 1,850 kg m −3 . Two kinds of clay mineral: kaolinite and montmorillonite, were used in the experiments. The sand is a uniform fine white quartz of a mean diameter around 0.23 mm and does not contain particles smaller than 63 µm. The mud (clay mineral) content varied from 0 to 14.9% for sandkaolinite mixtures and from 0 to 28% for sand-montmorillonite mixtures. The bed shear stress was estimated by the slope of the energy line which was calculated from the water surface profile. The critical shear stress was defined as the average of the bed shear stresses before and after the onset of erosion. The onset of erosion was reached when sediment is falling into the sediment trap and/or the water samples contain suspended sediment. As the formulae of Van Ledden (2003)

Experimental Data of Sharif (2003)
Sharif (2003) conducted a series of erosion tests of sand-kaolinite mixtures using a straight-through recirculating flume with a rectangular cross-section. The sand used in the experiments has a narrow particle size range from 180 to 212 µm. The median diameter of the sand is about 200 µm. The kaolinite used in the experiments is the industrial-grade pure kaolinite mineral. The mixtures were prepared by mixing the well-sorted fine sand, kaolinite and water with an electrical blender. The mud content varied from 0 to 100%. The mixtures were allowed to consolidate for 48 h before the experiments. For each mixture (or pure mud) bed, the critical shear stresses of the sediment at different depths from the bed surface were measured. Particles from all over the bed starting to erode from the bed was defined as the threshold condition for critical shear stress. The bed shear stress was estimated from the log-law for the vertical velocity distribution and a measured velocity at the center of the bed. Figure 5 shows the measured and predicted critical shear stress of sand-mud mixtures vary with mud content. Each subfigure represents the results of sediments of different mud contents for the same soil depth. For the formula of Van Ledden (2003), i.e., Eq. (1), p mcr is set as 15% and β v is set as 1.0. For the formula of Ahmad et al. (2011), i.e., Eq. (2), ζ is set as 0.15. These coefficients are given following the recommendations of Van Ledden (2003) and Ahmad et al. (2011) and are also used in other cases in this study. The best-fit value of A for Eq. (23) is 3.62 × 10 −6 J m −2 in this case. As shown in Figure 5, the calculated critical shear stress by the formulae of Van Ledden (2003)  There is a circular cut-off at the bottom of the duct so that a sediment container could be installed where a sediment core can be placed and pushed upwards. During the erosion testing, the surface of the sediment sample was horizontally and vertically leveled with the bottom of the duct. The bottom of the duct was covered with sandpaper of a roughness comparable to the applied sand fraction to decrease differences in roughness with the sample.
The sand-mud mixtures were prepared by mixing fine sand, silt and clay in different proportions. The median particle  diameters of the sand and silt used in the experiments are 180 and 30 µm, respectively. Two different clay minerals were applied: kaolinite and bentonite. A specific procedure was followed to generate reproducible and homogeneously mixed mixtures. First, sand, silt and clay fractions were oven-dried to disaggregate the material. Next, they were manually mixed for 10 min and placed in a cylindrical container with a removable bottom-lid. Small holes in the bottom and top-lid allow the passage of water and gas. In order to ensure the mixtures being 100% saturated, the containers with dry mixtures were placed in an exsiccator to remove air by lowering the pressure. Then the exsiccator was filled with CO2, after which the pressure was lowered again to replace enclosed air with CO2. Subsequently, the mixtures were left for 24 h in the exsiccator, in which a layer of water was present. The combination of the low pressure (reduced surface tension), 100% humidity and the attractive forces of the negatively charged clay particles enables water molecules to activate the clay particles. A layer of 10 cm de-aired and demineralized water was placed on top of the samples. Due to the difference between the atmospheric and the reduced pressure within the exsiccator, water percolated through the mixture thereby completing the saturation procedure.
To prevent consolidation, erosion tests were executed as soon as possible after the generation of the sample. A unidirectional flow generated by a recirculating pump was accelerated step by step (average duration of a step approximately 150-200 s) until the sample was eroded by a few mm. The volume of eroded sand was monitored at a sand trap downstream of the sediment sample, from which the erosion rate of the sand fraction could be derived. The erosion rate of the fine fraction is determined by dividing the time derivative of the continuously recorded turbidity by the surface area of the sample. The erosion threshold of a mixture was selected as the average abscissa of the extrapolated erosion rates assuming a linear relationship between the erosion rate and bed shear stress for both coarse and fine fractions. The bed shear stress was assumed proportional to the square of discharge in the flume, and the drag coefficient has been fitted so that the initiation of movement of the sand used in the experiments is consistent with the critical mobility parameter given by the Shields diagram.
As the mixtures were composed of sand, silt and clay and the clay contents are under 20% for all the mixtures, only Eq. (23) and the formula of Chen et al. (2018) are applied in this case. Figure 6 shows the measured and calculated critical shear stresses of mixtures varying with clay content. The best-fit value of A for Eq. (23) equals 1.16 × 10 −5 J m −2 for sand-silt-kaolinite mixtures and 5.93 × 10 −6 J m −2 for sand-silt-bentonite mixtures. As shown in Figure 6, the critical shear stresses calculated by Eq. (23) agree well with the measured values (Figures 6A,B), and the formula of Chen et al. (2018) tends to overestimate the critical shear stress at high clay contents (Figures 6C,D). Smith et al. (2015) conducted erosion tests of sand-mud mixtures using the Sedflume, which is a small-scale enclosed duct-type flume (similar to the Erodimetre introduced above). The sandmud mixtures were prepared by mixing varying fractions of mud with well-sorted quartz sand. The diameter of the sand ranges from 0.25 to 0.50 mm with the median diameter being 0.353 mm. Three different muds were used: kaolinite, mixed 80% kaolinite and 20% bentonite and a natural mud from the lower Mississippi River. The median diameters of the kaolinite, the kaolinite/bentonite mixture clay and the natural mud are 0.0038, 0.005, and 0.0127 mm, respectively. The natural mud has 3.5% organic content according to the test of loss-onignition.

Experimental Data of Smith et al. (2015)
Prior to mixing with sand, each mud was slurried with freshwater, and the slurry density was maintained at a density of about 1,400 kg m −3 . Then an appropriate mass of mud slurry was added to sand to achieve the targeted mud mass fraction. The mixture was then hand-mixed for 3-10 min until mud and sand were evenly distributed throughout the sample. The mixture samples were formed into 10-cm diameters cores by two methods. For fluid mixtures, samples were extruded from a bag into the bottom of a core to minimize gas entrapment. For plastic and granular mixtures, water was added to the mixture to achieve saturation, and the samples were lightly tamped into the core tube in approximately 0.5-1 cm lifts. The prepared sediment cores were allowed to consolidate for 30 days in a 4 • C cooler prior to the erosion test.
For each mixture sample, the erosion tests were conducted at five different shear stresses. The critical shear stress was defined as the bed shear stress corresponding to an erosion rate of 10 −4 cm s −1 and obtained by fitting a power law for the erosion rate. The bed shear stress in the Sedflume was determined by the average flow velocity in the flume according to the Darcy equation in which the friction factor was estimated by the Prandtl equation (i.e., both the bed of the flume and the surface of the sediment sample were assumed hydraulically smooth). Figure 7 shows the measured and calculated critical shear stresses varying with mud content. The formula of Chen et al. (2018) is not applied in this case since the stable dry bulk density of mud cannot be obtained based on the known data. The best-fit value of A for Eq. (23) equals 1.59 × 10 −6 J m −2 for sand-kaolinite mixtures, 6.42 × 10 −6 J m −2 for sand-kaolinitebentonite mixtures and 1.04 × 10 −5 J m −2 for sand-natural-mud mixtures. As shown in Figure 7, the formulae of Van Ledden (2003) and Ahmad et al. (2011) cannot capture well the varying behavior of the critical shear stress with mud content. Eq. (23) and the formula of Wu et al. (2017) reproduce well the critical shear stress of sand-kaolinite-bentonite mixtures ( Figure 7B) and sand-natural-mud mixtures ( Figure 7C). However, the latter two formulae can only reproduce the critical shear stress of sand-kaolinite mixtures of mud content lower than 15% and fail to capture well the varying trend of the critical shear stress of sand-kaolinite mixtures of high mud contents ( Figure 7A). The reason for the latter two formulae failing in capturing the varying trend of the critical shear stress of sand-kaolinite mixtures of high mud contents is probably because the nominal mud contents would be lower than the actual mud contents. For sand-kaolinite mixtures of high mud contents (i.e., fluid mixture samples in the preparing process), the uneven settling of the sand fraction and the mud fraction would occur during the preparing process of the samples. This would lead to that more mud is present in the upper layer of the sample. A modification is made here by assuming the nominal mud contents being 6% lower than the actual mud contents for the sand-kaolinite mixtures of mud content higher than 15%. The calculated critical shear stresses by Eq. (23) and the formula of Wu et al. (2017) after modification are also plotted in Figure 7A, which shows better agreements with measurements than before. The uneven settling of the sand fraction and the mud fraction is not significant in sand-kaolinite-bentonite mixtures and sand-natural-mud mixtures. This could be attributed to the mud fractions has much stronger cohesion for sand-kaolinitebentonite mixtures and sand-natural-mud mixtures, which could also be seen from the higher values of A.

VALIDITY OF THE DEVELOPED FORMULA IN PURE MUD
Nine sets of experimental data of pure mud were collected from literature to test the validity of the developed formula in pure mud. These muds include one group of pure clay mineral, two groups of estuarine mud and six groups of coastal mud. In each experiment, the critical shear stresses of the muds of different bulk densities were measured. The synopsis of the collected data and the data sources are listed in Table 2.
Equation (23) is then applied to the nine datasets. For each dataset, the best-fit value of A is used, which was obtained by the regression analysis of the data. The formula of Chen et al. (2018) and Eq. (4) (which was developed by Wu et al. (2017) for calculating the critical shear stress of pure mud) are also applied to each dataset if they are applicable. Figure 8 shows the application results of those formulae. As shown in Figure 8, the formula of Chen et al. (2018) and Eq. (4) only   reproduce well the critical shear stress of pure mud in a few cases. Eq. (23) gives reasonable estimates of the critical shear stress in all the cases, demonstrating the capacity of the developed formula in pure mud.

DISCUSSION
The Coefficient A The value of A is related to both the cohesion strength of the cohesive component in the sediment and the roughness of the bed surface. The cohesion strength is denoted by the Hamaker constant, which is a function of the particle itself including the mineral composition, the particle shape and the surface roughness, and the pore water environment including the pH value, the sort and concentration of ions. As the mineral composition and the pore water environment of sediment usually vary from site to site, the value of Hamaker constant and the value of A are supposed to be site-specific. The surface roughness of sand-mud mixtures depends on the network structure of mixtures, the aggregation of mud and the occurrence of erosion events (Mitchener and Torfs, 1996;Le Hir et al., 2008;Perkey et al., 2020). Different experimental results of surface roughness of sand-mud mixtures were often reported (Baas et al., 2013;Das et al., 2019). A widely accepted quantitative understanding of the roughness parameter of sand-mud mixtures and its effect on the drag and lift coefficients have not been available. Therefore, A is treated as an empirical coefficient in this study with its value being determined by the measured erosion thresholds. The threshold of erosion in most existing tests of sandmud mixtures and pure mud was usually determined by visual observation which is subjective and highly empirical. A universal criterion for the erosion threshold of sediment, especially of cohesive sediments, is not available currently. Therefore, the different criteria for determining the threshold of erosion would also have an effect on the value of A. This makes it more challenging to develop an expression for A.
The best-fit values of A for each dataset of sand-mud mixtures and pure mud are presented in Table 3. As shown in Table 3, the value of A is generally on the order of magnitude of 10 −6 or 10 −5 J m −2 . Specifically, the value of A is the range of 1.59 × 10 −6 −1.16 × 10 −5 J m −2 for kaolinite and sand-kaolinite mixtures; 5.93 × 10 −6 -9.58 × 10 −6 J m −2 for sand-bentonite mixtures; 4.89 × 10 −6 J m −2 for mixtures of sand and kaolinite and bentonite mixed clay; 2.86 × 10 −6 -1.04 × 10 −5 J m −2 for pure mud and sand-mud mixtures.
Although the value of A is site-or sediment-specific, a general or practical value is expected, with which Eq. (23) could give a reasonable estimate of critical shear stress of most of the common sediments. Based on the collected data of sand-mud mixtures and pure mud, a general value of 3.97 × 10 −6 J m −2 is recommended, which was obtained by the following method. The value of A was allowed to increase from 1.0 × 10 −7 J m −2 to 1.0 × 10 −4 J m −2 with a step of 1.0 × 10 −8 J m −2 . For each value of A, Eq. (23) was applied to all the collected datasets, and the logarithmic rootmean-square error of the calculated critical shear stresses was computed. The logarithmic root-mean-square error is used as a statistical indicator to evaluate the performance of the regression model for the quantity of the measured critical shear stress varies in several orders of magnitude. The general value of A is obtained when the logarithmic root-mean-square error of the calculated critical shear stresses reaches its minimum value. The logarithmic root-mean-square error, denoted as log E rms , is defined by: where τ cr,c,i and τ cr,m,i are the calculated and measured critical shear stresses, respectively; i is the data index and M is the total number of the collected data. The comparisons of Eq. (23) and the measured critical shear stresses in all the datasets are shown in Figures 9A,B. For Figure 9A, the best-fit value of A (see Table 3) is used for each dataset; and for Figure 9B, the general value of A, 3.97 × 10 −6 J m −2 , is used for all the datasets. The comparisons of the formulae of Wu et al. (2017) and Chen et al. (2018) with measurements are also shown in Figures 9C,D. As shown in Figure 9, Eq. (23) with best-fit values of A has the best prediction performance. Eq. (23) with the general value of A also has a good prediction performance, better than the formulae of Wu et al. (2017) and Chen et al. (2018). The average relative error (E), defined as E = τ cr,c,i /τ cr,m,i − 1 /M × 100% , is 22.9% for Eq. (23) with best-fit values of A, 34.3% for Eq. (23) with the general value of A, 38.3% for the formula of Chen et al. (2018) and

Effects of the Network Structure of Sand-Mud Mixtures
Equation (23) indicates that the dry bulk density of the mud component is a more direct factor affecting the erodibility of a mixture than the dry bulk density of the mixture. The variation of the critical shear stress of sand-mud mixtures with mud content is mainly caused by the varying dry bulk density of the mud component in the mixture. These findings are consistent with previous observations by Migniot (1989); Waeles et al. (2008), Dickhudt et al. (2011); Le ), Wu et al. (2017, and Chen et al. (2018). As the network structure of sand-mud mixtures affects the consolidation of sand-mud mixtures, the network structure of sand-mud mixtures would have effects on the dry bulk density of the mud component and also the mixture critical shear stress. Here the effects of the network structure of sand-mud mixtures are investigated using the experimental data of Sharif (2003). Figure 10 shows the dry bulk density of the mud component ( Figure 10A) and the critical shear stress calculated by the developed formula ( Figure 10B) varying with mud content in the case of Sharif (2003). A is taken the value of 3.62 × 10 −6 J m −2 in the calculations. All the sediments in the experiments were allowed to consolidate for 48 h before the erosion test. The dry bulk density and the critical shear stress of sediment at different depths from the bed surface were measured. Each solid line with markers in Figure 10 represents the results of sediment samples of different mud contents at the same depths. Figure 10A shows that the dry bulk density of the mud component increases markedly with increasing mud content from 0 to 10%. With further increase in mud content, the dry bulk density of the mud component first continues to increase and reaches a maximum at an optimum mud content, after that decreases gradually in a slight rate until 100% mud content. The optimum mud content seems to be in the range of 30-50% in this dataset. The dips of the dry bulk density at 40% mud content may be caused by uncertainty in measurements. As generally predicted by the developed formula, the critical shear stress varies with mud content in a similar way as the dry bulk density of the mud component (Figures 10A,B).
As shown in Figure 10A, the dry bulk density of the mud component does not change with depth for the mud content below 10% and varies with depth when the mud content equals or exceeds 20%. This indicates that the critical mud content for being able to form a sand skeleton is between 10 and 20%. Because when a sand skeleton is formed, the mud particles are in the voids of the sand particles. For such a condition, the mud particles are protected by the sand skeleton, and the sediment depth will have no effect on the consolidation of the mud component. The mud component only consolidates under its self-weight, and therefore, the more mud fraction, the larger the dry bulk density of the mud component. When the mud content is beyond the critical mud content, the complete sand skeleton cannot be formed. For this case, the mud component consolidates under the influence of the effective weight of the overlying layer. Therefore, for the same consolidation time, the deeper sediment depth, the higher the degree of the consolidation of the mud component.
The behavior of the dry bulk density of the mud component increasing first and then decreasing gradually when the mud content is beyond the critical mud content is not wellunderstood. A reasonable explanation is proposed here. When the mud fraction is not very high, although a complete sand skeleton cannot be formed, some sand particles could still contact each other and form an incomplete skeleton. For this condition, the mud component and the sand component share the load from the effective weight of the overlying layer of sediment. The mud component consolidates under the received load and the self-weight. The more mud fraction means less help from the sand particles and more load received by the mud component. Therefore, a higher mud content will lead to a faster consolidation of the mud component in a mixture. However, if the mud fraction is high enough that all the sand particles lose contact and are dispersed in the mud matrix. For this condition, the mud component consolidates under the self-weight and the load from the overlying layer. As the sand could assist drainage and the density of the sand particles is higher than that of mud blocks, the increase in the mud content will not benefit the consolidation of the mud component. The more mud fraction and the less sand fraction lead to a lower consolidation rate. Under these two opposing influences, the dry bulk density of the mud component in sand-mud mixtures of the same consolidation time will first increase then decrease.

CONCLUSION
The erosion threshold of sand-mud mixtures has been studied by a theoretical analysis of the momentum balance of a sand particle or a mud parcel under the initial motion condition. A formula for the critical shear stress of sand-mud mixtures was developed, which also applies for pure sand and mud. The developed formula is expressed as a function of the primary particle diameters of the sand and mud components, mud content and the dry bulk density of sediment. Compared with the existing formulae for sand-mud mixtures, the developed formula is much easier for application as it has few coefficients and does not involve the critical shear stress of pure mud and the stable dry bulk density of mud in its expression.
The developed formula can be simplified to the widely used formula for pure sand when the mud content is 0 and has been successfully tested by four sets of experimental data of sand-mud mixtures and nine sets of experimental data of pure mud, which demonstrates the capability of the developed formula in pure sand, pure mud and sand-mud mixtures. The value of A in the developed formula is supposed to be site-specific. According to the collected experimental data, it is generally on the order of magnitude of 10 −6 or 10 −5 J m −2 . A value of A, 3.97 × 10 −6 J m −2 , is recommended for general sand-mud mixtures and mud. But it is emphasized that the value of A could be optimized when the developed formula is applied to a specific site.
The developed formula shows that the erosion threshold of sand-mud mixtures is highly affected by the dry bulk density of the mud component. The varying dry bulk density of the mud component in a sand-mud mixture is the main reason for the variation of its critical shear stress with mud content. By analyzing the consolidation behavior of the mud component in mixtures, the effects of the network structure of sand-mud mixtures on the dry bulk density of the mud component and the critical shear stress of sand-mud mixtures were investigated.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
DC and CZ contributed to the conception of the study, the theory development, and manuscript preparation. JZ and YW contributed significantly to the theoretical analysis. DG helped perform the analysis with constructive discussions. YL contributed to the analysis of the data for the work and helped with revising the work. All authors contributed to the article and approved the submitted version.