Detachment of dilatant soil due to high hydraulic shear stresses explained

When soil surfaces are exposed to high shear stresses induced by high-velocity water flows, they erode. Erosion consists of the detachment, entrainment and transport of soil particles. When exposed to high shear stresses the process of soil detachment is no longer given by the pick-up of individual particles but by the continuous shear failure of layers of soil. The understanding of how soil properties influence the detachment rate of soil is lacking. This paper presents how the bed shear stress, initial and critical porosity, hydraulic conductivity, and density influence the detachment rate of dilatant soils which are subjected to high-velocity flows. Based on the constitutive mass and momentum balance equations it is hypothesized that the exchange of mass and momentum during soil detachment corresponds with that situation for which the dilatancy induced shear resistance of the soil is maximum. The hypothesis is validated against high-velocity erosion experiments on sand.


Introduction
Shear stresses induced by high-velocity water flows past soil surfaces, give erosion. Erosion consists of the detachment, entrainment and transport of particles. The development of accurate and widely applicable predictive relations for quantifying erosion rates is hindered by a lack of understanding of the process of soil detachment. Hanson and Hunt (2007) and Morris, Hanson, and Hassan (2008) noted that the detachment rate or pick-up flux is dependent on material texture, compaction moisture content, and compaction energy. When dilatant soils are subjected to high surface shear stresses, soil detachment is no longer characterized by the pick-up of individual particles but by the continuous shear failure of layers of soil (Van Rhee, 2010). During soil detachment, particles and water are exchanged at the bed. With the exchange of particles and water, momentum is also exchanged. Above the depth at which the shear stresses equal the shear strength, the bed begins to move (Iverson, 2000;Iverson, Reid, & LaHusen, 1997;Takahashi, 1981Takahashi, , 2009. For layers to shear, soil needs to dilate. The corresponding increase in pore volume requires an inflow of water into the bed, which in turn requires a pressure gradient to be present. As the pressure on the bed is hydrostatic, pore water pressures in the soil need to decrease for a water flow to initiate. As the effective stresses in the soil consist of the soil pressures minus the water pressures, a reduction in pore water pressures gives an increase in effective stresses and shear resistance (Van der Schrieck, 2006). A low initial porosity requires a high degree of dilation, which results in a high shear resistance. A low permeability requires a large pressure gradient to accommodate the inflow, and consequently also results in a high shear resistance. When a soil surface is exposed to high flow velocities, the rate of detachment therefore depends on the initial porosity with respect to the critical porosity, and the hydraulic conductivity of soil (Bisschop, Miedema, Visser, Keetels, & Van Rhee, 2016;Van Rhee, 2010). The critical porosity is thereby the porosity for which shear is possible without a change in pore volume. For sand, the influence of dilation already becomes significant for shear stresses corresponding with flow velocities in excess of approximately 1 m s −1 (Bisschop et al., 2016). The effects of dilation correspond well with observations that soil erodibility depends on material texture, compaction moisture content, and compaction energy (Hanson & Hunt, 2007;Morris et al., 2008). Although the impact of dilation has been identified (Van Rhee, 2010;Winterwerp, Bakker, Mastbergen, & van Rossum, 1992), accounting for the effects of dilation has been challenging. Mastbergen and Van den Berg (2003); Van Rhee (2010); Winterwerp et al. (1992) published sediment transport formulas for flow velocities up to 3 m s −1 . These models give better predictions for the pick-up flux than Van Rijn's pick-up function (Bisschop et al., 2016;Van Rijn, 1984) but still have a predominant empirical nature. A disadvantage of using empirical equations is that ample data is required to validate the equation and that the range of applicability is limited by the test conditions.
Here it is shown that the rate of detachment of dilatant soils corresponds with a maximum in shear resistance induced by dilatancy, which is a function of the porosity, hydraulic conductivity, and soil density. Section 2 of this paper outlines the mass and momentum balance equations which describe the effect of dilation on the exchange of water and particles at the bed. The method for predicting detachment rates is validated against experimental data in Section 3. A discussion of the results and the conclusions are given in the Sections 4 and 5.

Methodology
This section deals with the behaviour of a bed of dilative sediment which is subjected to shear stresses due to fast flowing water. Shear stresses applied to a bed are transferred to lower lying soil layers. Particles, close to the bed, are mobilized when the shear stresses exceed the shear resistance (Takahashi, 2009). For a horizontal bed, the bed shear strength is defined by the Coulomb's friction which is given by: where τ c,zx denotes the shear strength component in the xdirection, σ zz denotes the effective stress component in the z-direction, φ is the angle of internal friction, and c is the cohesion of the bed. Dilation of the soil is required for particles and soil layers to move over each other (Fig. 1). The associated increase in pore volume requires an inflow of water into the bed for which a pressure gradient must be present over the bed.
The effective stress components in soil are given by the total stress minus the pore water pressure P. For a constant total stress, a reduction in pore water pressures gives an increase in effective stress and shear resistance (Eq. 1). The process of detachment of dilatant material is therefore a balance between the shear stress and the increase in effective stresses and shear resistance induced by dilatancy. At the top of the bed the exchange of water and particles indicates an exchange of momentum. When the rate of momentum exchange between the bed and the particles equals the applied shear stress, a balance is achieved (Iverson et al., 2011).
The following sections discuss the constitutive equations that describe this balance in x, z, t coordinates. The x-coordinate has thereby been taken parallel to the bed-water interface. The z-coordinate direction has been defined in the direction perpendicular to the bed-water interface (Fig. 2). t denotes the time coordinate. The positive vector components for the velocity u and w respectively, coincide with the positive x-and z-coordinate directions.

Mass balance equation
Assuming incompressible particles and incompressible water the mass balance equation for water in a porous system is given by: where u w denotes the horizontal flow velocity component, w w denotes the vertical flow velocity component of the water in Similarly the mass balance equation for solid particles in a porous medium is given by: where u p denotes the horizontal flow velocity component and w p denotes the vertical flow velocity component of particles in the porous medium. Adding Eqs (2) and (3) gives: The last two terms in Eq. (4) describe a change in porosity. This has been accounted for by replacing these terms by ∂ /∂t. Here denotes the volume strain, for which (1/(1 − n))(∂n/∂t) is substituted (Verruijt, 2006). The terms n(u w − u p ) and n(w w − w p ), respectively, denote the specific discharge components q x and q z , which are related to the pressure gradient by Darcy's law. Substituting the expressions for the specific discharges in Eq. (4) gives: The bed location continuously changes as soil layers are sheared off. To describe the continuous displacement of the bed a moving ξ , η coordinate system is defined whereby ξ = x and η = z − vt (Jensen & Finlayson, 1980;Verruijt, 2006). Here v (m s −1 ) denotes the rate of displacement of the bed in the z-coordinate direction. The thickness of the sheared layers is small compared to the characteristic length scale of the bed surface. This allows for a one-dimensional description of the erosion problem. The corresponding 1D mass balance equation for water (Eq. 2) now becomes: The depth over which shear stresses exceed the shear strength, denoted by Coulomb's friction, is from this point onwards denoted by d.
During dilation, the flow of water towards the shear zone induces a pressure drop. Deep in the bed, pore pressures are likely to be unaffected. At the shear plane, pressures are minimum. On both sides of the shear plane consequently a pressure gradient is present, which is 0 at the shear plane. The flow rate coming from the soil below the shear plane is limited by the degree of contraction, which is assumed negligible over the short time duration over which thin soil layers shear off. Consequently, water is assumed to predominantly enter the bed through the bed surface. In line with this the velocity component w w is assumed 0 at depth d. Integrating the mass balance equation for water over η from η = 0 to η = d gives for the mass balance of water: Rewriting Eq. (3) in terms of a moving coordinate system and integrating from η = 0 to η = d along the same lines as done for the mass balance equation of water (Eqs 6 and 7) gives for the mass balance of particles: whereby it has been assumed that w p (1 − n)| η=d = 0. The mass balance of the total porous medium follows from rewriting Eq. (5) in terms of a moving coordinate system, and integrating from η = 0 to η = d. This then leads to: Here n| η=d is assumed to equal the initial porosity n 0 , and n| η=0 is assumed to equal the porosity of the fully dilated material n loose (Fig. 3) (Van der Schrieck, 2006).

Momentum balance equations
The general momentum balance equations for soil consist of the momentum balance of particles and that of pore water. Below the momentum balance equations for an infinitely long slope are given whereby changes in the x-coordinate direction are 0. The direction of dilation is thereby in the negative z-coordinate direction. The momentum balance equation in the x-coordinate direction for erosion of a surface to which a uniform shear stress is applied, is given by: where ρ p , ρ w respectively, denote the density of soil particles and the density of water, g is the gravitational constant, K s is the hydraulic conductivity of the soil, and τ zx is the horizontal shear stress component. The surface shear stress induces dilation of the soil, during which momentum is transferred via the pore water pressures to the soil skeleton. The momentum balance equation describing the exchange in momentum in the z-coordinate direction via the pore water pressures is given by: Figure 3 Sediment balance at an eroding surface, where v denotes the bed displacement rate, n loose denotes the porosity of the fully dilated material, and n 0 denotes the initial porosity of the bed. The time steps t and t + t denote the change in location of the bed level (after Van der Schrieck, 2006) where P denotes the pore water pressure. Currently the nonlinear advective acceleration terms are in a format which is not easily integrated. To facilitate integration, Eqs (7) and (8) have been substituted in Eq. (10) under the assumption that n is spatially constant over d (Van der Schrieck, 2006;Van Rhee, 2010;Verruijt, 2006). This assumption is in line with conclusions of Casagrande (1936) who observed that samples subjected to large deformations reach the same porosity (Van Rhee, 2010). Hence over the depth d the shear induced dilation of soil continues until the critical porosity is reached at which particles are able to shear over each other. Replacing the fixed x, z, t coordinate system with the moving ξ , η, t coordinate system finally gives: where τ ηξ denotes the horizontal shear stress component in the moving coordinate system. The same process of substituting the mass balance equations and transforming the fixed coordinate system into a moving coordinate system has been applied to the vertical momentum balance equation (Eq. 11), leading to: In line with the assumptions made to arrive at the integrable form of the momentum balance equations, dilatancy is assumed to cause an equal rate of increase in the pore spaces over the interval [0, d). At the depth d dilation is 0 and u w = w w = u p = w p = 0. Over the shear depth it is assumed that u w is equal to u p . Integrating Eqs (12) and (13) under these conditions from η = 0 to η = d gives: and: The mass balance described by Eqs (7) and (8) state that at the surface (η = 0), w w | η=0 = v((n loose − n 0 )/n loose ) and w p | η=0 = −v((n loose − n 0 )/(1 − n loose )). Substituting these relations in the momentum balance equations gives the following integrated momentum balance equation in the ξ -coordinate direction expressed in terms of rate of displacement of the boundary v: for the momentum balance equation describing the exchange of momentum in the η coordinate direction this gives: The last term in Eq. (17) denotes the boundary condition for the pressure on top of the soil given by P| η=0 in Eq. (15). D| η=0 thereby denotes the water depth on the soil and γ w = ρ w g denotes the specific weight of water. The first term on the second row of Eq. (17) equals the total friction experienced by the flow over the failure depth d. The average infiltration velocity over d thereby equals half the infiltration velocity at the surface in order for degree of dilatancy to be spatially constant. The rate of displacement of the soil body is of the order of mm s −1 , and the failure depth is of the order of mm. With K s being of the order of 1 × 10 −3 m s −1 or lower, the effects of the friction far outweigh the effects of the advective acceleration indicated by the terms containing v 2 . The effective stresses follow from the soil pressure minus the water pressure P| η=d , which follows from Eq. (17). The soil pressure at the failure depth is assumed to equal σ ηη = γ s d + γ w D| η=0 . Subtracting the water pressure leaves the effective stress component σ ηη | η=d , which is given by: Here γ s denotes the specific weight of the soil, which follows from γ s = n 0 γ w + (1 − n 0 γ p ), where γ p is the specific weight of the particles. The failure depth is given by that depth at which the absolute value of the shear stress component in the ξ direction equals the absolute value of the shear strength component in the ξ direction (Verruijt, 2001 Equation (19) relates the failure depth d to the rate of displacement of the boundary v. As v increases, d decreases.
The horizontal near-bed flow velocity component u w = u p has thereby been approximated by the shear velocity (Schiereck, 2004), which is given by: In Eq. (19) both the failure depth d and displacement rate v are unknown. An extra relation is hence required to close the equations. This extra relation, on which is hypothesized below, is the main contribution of this paper. Erosion relations often relate the erosion rate to an excess shear stress. In line with this it has been evaluated how the average shear resistance varies with the erosion rate v for a constant bed shear stress τ ηξ | η=0 . The degree of dilation over d is assumed constant, indicating that w w and w p decrease linearly from η = 0 to η = d as given by: The average pore water pressure P, induced by dilation, follows from twice integrating with respect to η and dividing by d.
Applying the pressure boundary conditions that the normal pressure gradient dP/dη| η=d = 0 and that the pressure P| η=0 = 0 then leads to: The average pore water pressure over the failure depth follows from substituting η = d. The average effective stress component σ ηη follows from the total stress minus the average pore water pressure. The overall shear resistance of the soil is maximum when the averaged effective stress is maximum. The expression for the average effective stress component in ξ direction over d is given by: where d follows from rewriting Eq. (19) into: Substituting Eq. (24) in Eq. (23) gives a relationship for the average effective stress as a function of the detachment rate of the soil surface v. This is thereby a key step in the derivation of the process based pick-up relation. Figure 4 shows this relationship for a constant shear stress, constant porosities n 0 and n loose and three different values for the hydraulic conductivity. Varying any of the other parameters in Eqs (24) and (23) gives similarly shaped curves. Figure 4 forms the basis for the main hypothesis of this paper, namely that during erosion, the bed displacement rate v corresponds to a maximum in bed shear resistance induced by dilatancy. This maximum corresponds to a minimum in excess shear stress, or momentum, available for the exchange of particles and water at the bed.

Results
Here the main hypothesis of the paper is validated against the dataset published by Bisschop (2018), who performed erosion experiments under high flow velocities on different sand types. Bisschop (2018) performed the experiments using a closed circuit of pipes through which a sand-water mixture could be pumped (Fig. 5). Within the circuit a parallel section allowed for the flow to bypass the measurement section in which a sand bed was prepared. Prior to the erosion tests a sand-water mixture of the required density was pumped around via the bypass. Once the required flow velocities were obtained, the flow was redirected through the measurement section to initiate the erosion process. The pump circuit was equipped with the following instruments to monitor the erosion process: • Radioactive density meter to monitor the mixture concentration • Conductivity probes in the measurement section to monitor the change in bed level during erosion • Pressure gauges to determine the pressure gradients over the measurement section, and the density of the sand-water mixture. • Electro-magnetic flow meters to determine the total discharge • A temperature sensor The properties of three types of dilatant sand beds tested by Bisschop (2018) are given in Table 1. Here n min and n max , respectively, refer to the minimum and maximum porosity. For validation, n loose has been assumed equal to n max . The parameters d 10 , d 50 , and d 90 , respectively, refer to the diameters at which 10%, 50%, or 90% of the sample's mass is comprised of particles with a smaller diameter. The uniformity of the particle distribution is indicated by the uniformity coefficient c u .
Between tests, the initial porosity, flow density, and flow velocity were altered to investigate the impact of these parameters on erosion. Initial porosities (n 0 ) thereby ranged from 0.356 to 0.488, and flow velocities ranged from 1.3 to 6.3 m s −1 . Mixture densities of the flow over the bed were varied between 1030 and 1390 kg m −3 to determine their effect on erosion. The mixture densities were measured using the radioactive density meters. The bed shear stress was provided by Bisschop (2018) for each of the tests. Bisschop (2018) measured the density of the bed using conductivity probes. Erosion rates were estimated from the drop in density of the sand bed just prior to the onset of erosion. During tests, both pick-up and deposition of the sediment occurred at the same time. Bisschop (2018) converted the measured erosion rates into a pick-up flux by accounting for the rate of deposition. Predictions for the detachment rate v have therefore been converted into a pick-up flux E via E = v(n 0 ρ p ) (Bisschop, 2018). The match between the predicted pick-up flux and measured pick-up flux for the different sand types are given in Figs 6-8. The spread in predictions shown in Figs 6-8 originates from varying the input parameters like flow density, and the initial and critical porosity. The difference between the highest and lowest predicted pick-up fluxes for comparable shear stresses is smaller than the difference in measured pick-up fluxes. During the experiments the bed shear stress was obtained from the pressure drop over the test section (Bisschop, 2018). Difficulties with predicting the rate of change in flow height between two time steps and uncertainties regarding the effects of the pick-up flux on variations in flow density were found to induce relative errors of τ ηx| η=0 /τ ηξ | η=0 in the range of 0.05 and 0.5 (Bisschop, 2018). This could explain the  Research Vol. 59, No. 1 (2021) Detachment of dilatant soil due to high hydraulic shear stresses explained 57  Bisschop, 2018) larger variation in measured pick-up fluxes than in the predicted pick-up fluxes for comparable shear stresses. Bisschop (2018) performed the erosion experiments for various sediment concentrations. Figure 9 shows how changes in mixture density of the flow affected the bed shear stress during the experiments. Figures 10-12 show the correspondence between the measurements and predictions for how the pickup flux varied with mixture density for the various sand types. These figures show that the increase in shear stress with mixture density is higher than the increase in pick-up flux with mixture density. Morris et al. (2008) noted that in order to improve the accuracy of erosion predictions, the understanding of the role  Bisschop, 2018) of material type, compaction effort, and moisture content on soil detachment needs to be improved. The initial porosity depends on the compaction effort and compaction moisture content. The inflow rate of water into the bed, which has a stabilizing effect, is a function of the initial moisture content of the bed. Finally, the material type influences the hydraulic conductivity. Consequently, the role of material type, compaction efforts, and moisture content influences can be explained from the effects of dilation on erosion. From Eqs (23) and (24) follows that a densely packed soil erodes slower than a loosely packed soil. A low permeable soil thereby also erodes slower than a more permeable soil (Fig. 4) (Bisschop, 2018). Figures 6-8 show that the predictions in pickup flux that correspond with a maximum in shear resistance induced by dilatancy, are in good agreement with the measurements. The smaller spread in measurements than predictions can thereby be attributed to measurement uncertainties identified by Bisschop (2018). Figures 9-12 illustrate the effects of the density of the flow on the erosion rate. The observed increase in erosion rate with density, as shown in Figs 10-12, is caused by the increase in bed shear stress with the density of the flow. Besides inducing a higher bed shear stress, the density of the flow may also influence the momentum exchange at the bed, denoted by the acceleration terms in the depth integrated momentum balance equation (Eq. (19)). Substituting the density of mixture flow over the bed for the density of water in the acceleration terms of Eq. (19), however, had a negligible effect on the prediction of the erosion rate. The increase in pick-up flux with density is consequently expected to be solely attributable to the increase in bed shear stress. Dilatancy has been shown to have a significant impact on the erosion rates (Van Rhee, 2010). Erosion experiments under high flow velocities, for various bed level and water level gradients, are recommended to expand the applicability of the theoretical analysis.

Discussion
Underlying the mass and momentum balance equations are the assumptions that the degree of dilation is constant above the depth at which the shear stresses equal the shear strength, and that water only enters the soil through the bed surface. The first of these assumptions was required to arrive at an integrable form of the momentum balance equations and is in line with the conclusion of Casagrande (1936), who observed that samples subjected to large deformations reach the same porosity. The assumption that water enters the soil through the bed was necessary due to a lack of information on the flow of water from the soil towards the shear zone. The agreement between the measurements and predictions indicates that these assumptions are acceptable. An extension of the method to include the effects of gradients in bed level or water level is recommended to extend the range of applicability of the method.

Conclusions
This paper explains how for dilatant soils subjected to high bed shear stresses, the pick-up rate is related to the process of dilation. As illustrated by Fig. 4, a maximum in shear resistance, induced by dilatancy, is found for a specific sediment pick-up rate. This pick-up rate corresponds with a minimum in excess shear stress, or momentum, available for exchange of particles and water at the bed. The pick-up rate is consequently a function of the bed shear stress, initial and critical porosity, hydraulic conductivity, and density. Studies on the effects of gradients in bed slope and water level on erosion are recommended to extend the applicability of this relation.

Funding
This research was funded by Stichting voor de Technische Wetenschappen (NWO TTW) (Project Number 13861).