Modelling suspended sediment in environmental turbulent fluids
Abstract
Modelling sediment transport in environmental turbulent fluids is a challenge. This article develops a sound model of the lateral transport of suspended sediment in environmental fluid flows such as floods and tsunamis. The model is systematically derived from a 3D turbulence model based on the Smagorinski large eddy closure. Embedding the physical dynamics into a family of problems and analysing linear dynamics of the system, centre manifold theory indicates the existence of slow manifold parametrised by macroscale variables. Computer algebra then constructs the slow manifold in terms of fluid depth, depthaveraged lateral velocities, and suspended sediment concentration. The model includes the effects of sediment erosion, advection, dispersion, and also the interactions between the sediment and turbulent fluid flow. Vertical distributions of the velocity and concentration in steady flow agree with the established experimental data. Numerical simulations of the suspended sediment under large waves show that the developed model predicts physically reasonable phenomena.
Contents
1 Introduction
Environmental turbulent fluids, such as rivers, floods and tsunamis, always carry amounts of sediment. For example, Figure 1 shows the Yellow River in China which is famous for carrying large amounts of sediment. Modelling the sediment in these environmental fluids is important for studying and predicting changes of the morphology and topography. We aim to develop a model to appropriately model the suspended sediment in turbulent fluid flows via systematic resolution of the physical processes.
Our modelling, which is based on dynamical systems theory instead of conventional depthaveraging, resolves outofequilibrium interactions between the varying turbulence and suspended sediment. Most previous work studied suspended sediment in uniform flows (Hunt, 1954; van Rijn, 1984; Celik & Rodi, 1988; Fredsoe & Deigaard, 1992, e.g.) or by depthaveraging flow and sediment equations (Wu et al., 2000; Pittaluga & Seminara, 2003, e.g.). We explore the implications of changing the theoretical base from depthaveraging to a slow manifold of the turbulent Smagorinski large eddy closure.
Cao & Roberts (2012) initially developed a 2D lateral model for environmental fluid flows, derived from a 3D turbulence model based on Smagorinski large eddy closure. This model includes the effects and interactions of inertia, advection, bed drag, gravitational forcing and turbulent dissipation with minimal assumptions. The innovation here is that the turbulent modelling and dynamics includes and interacts with the suspended sediment transport Cao (2014). A slow manifold is found for the outofequilibrium dynamics of the coupled turbulent sediment system. We choose to parametrise the slow manifold in terms of emergent depthaveraged quantities. The evolution of these depthaveraged quantities on the slow manifold governs the dynamics of the suspended sediment in the turbulent fluid flows.
Consider a turbulent flow of depth flowing along a bed with a mean slope and carrying sediment. Sections 2.1–2.2 detail equations of the Reynoldsaveraged Navier–Stokes pdes, the advectiondiffusion equation and boundary conditions on the free surface and the mean bed. Section 2.3 uses the dynamical systems theory of centre manifolds (Roberts, 1988; Potzsche & Rasmussen, 2006, e.g.) to analyse the governing equations and derive the following nondimensional suspended sediment model of the horizontal evolution of the depth , the lateral depthaveraged velocities and and depthaveraged concentration :
(1a)  
(1b)  
(1c)  
(1d) 
where the mean local flow speed, constant is the falling velocity of the sediment, and constant is an equilibrium reference concentration on the mean bed . The effective momentum equations (1b)–(1c), and the more refined version (28b)–(28c), include the effects of gravitational forcing, bed drag, selfadvection, turbulent dissipation, and sediment induced flow. The sediment pde (1d) includes the sediment erosion and deposition, and the advection. Although this model is expressed in terms of depthaveraged lateral velocities and concentration, they are derived not by depthaveraging, but instead by systematically accounting for interaction between vertical profiles and horizontal gradients of the velocity, the concentration, the stress, bed drag, lateral space variations and bed topography.
That the systematically reduced model (1) has many terms so close to established models is a testament to the robustness of conservation of mass and momentum principles that underly traditional derivations. Additionally, our dynamical systems approach resolves finer microscale and outofequilibrium interactions that are active in more physically complicated systems.
Linear analysis in section 2.4 finds the spectrum supports the existence of a slow manifold. The computer algebra of Appendix A then constructs the slow manifold of the system. Then section 2.5 derives the reduced model (1) of the turbulent flow and suspended sediment on the slow manifold.
Section 3 discusses the predicted vertical distribution of the velocity and concentration fields in steady flow. Agreement with established experimental data (Schultz & Flack, 2007, 2013; Celik & Rodi, 1988, e.g.) indicates that our model is reasonable. Section 4 numerically explores the suspended sediment under large waves by the comprehensive suspended sediment model (28d) coupled with the turbulence model (28a)–(28c). The numerical results indicate that our model reasonably describes the dynamics of the suspended sediment in turbulent flows.
2 Construct the sediment transport model
This section describes the derivation of the reduced models for turbulent flow and suspended sediment. First, section 2.1 details the 3D continuity, Navier–Stokes and advectiondiffusion equations of the turbulent fluid flows and suspended sediment, whereas section 2.2 records the boundary conditions of the flow and the sediment on the free surface and the mean bed. Second, section 2.3 embeds these equations in a family of equations with modified tangential stresses on the free surface to establish the existence of an appropriate slow manifold. The linear analysis of section 2.4 supports the emergence of the slow manifold from the dynamics in the system. The computer algebra of Appendix A handles the details of the construction of the slow manifold model that is summarised in section 2.5.
2.1 The governing equations of the flow and sediment
Consider the turbulent flow flowing along a bed carrying sediment. This work only considers the suspended sediment and neglects the bed load transport on the mean bed. Figure 2 depicts the diagram of the suspended sediment in the turbulent flow. Define a coordinate system with for the lateral directions, and for the direction normal to the mean slope.
The fluid of depth flows down the sloped mean bed at the turbulent mean velocity ; the velocity vector in the directions, respectively. The term ‘mean bed’ refers to the smooth average bed over an ensemble of turbulence and bed roughness realisations: the physical bed in any one realisation is envisaged to be rough, just like the physical fluid velocity field would have rapid spatial variations in any one realisation. Denote the turbulent mean pressure field by . The suspended sediment has a turbulent mean concentration (volume fraction). The mean bed has an overall slope in the direction.
We assume the particles of the suspended sediment have small sizes, and all the particles of the suspended sediment have the same falling velocity, namely that of a sphere of diameter .
The nondimensional governing partial differential equations for the incompressible, three dimensional, turbulent mean fluid fields are the Reynoldsaveraged continuity and momentum equations,
(2)  
(3) 
and for the suspended sediment is the advectiondiffusion equation,
(4) 
where is the falling velocity of the sediment particles, is the eddy viscosity, and the vector is the direction of the nondimensional gravity.
The suspended sediment influences the fluid turbulence. We assume the mixing density of the fluid and suspended sediment satisfies
(5) 
where is the sediment density, the density of the fluid and is the relative density.
The pdes (2)–(5) are nondimensional and derived upon using the characteristic depth of the turbulent fluid as the length scale, the long wave speed as the velocity scale, and the fluid density as the reference density. Thus, the mixing density .
The variable is the turbulent mean deviatoric stress tensor, which is approximated using the eddy viscosity . We use the Smagorinski eddy closure to approximate the turbulence stresses. Cao & Roberts (2012) and Georgiev et al. (2009) detailed this eddy closure and expressed the mean deviatoric stress tensor
(6) 
where the magnitude of the second invariant of the strainrate tensor satisfies .
The falling velocity is related to the mean particle effective size , the relative density and the gravity . We set the falling velocity (Fredsoe & Deigaard, 1992, e.g.)
(7) 
The drag coefficient for the large grain Reynolds number of natural sands, typically (Fredsoe & Deigaard, 1992, e.g.). Typically, a quartz sediment has a relative density (Chanson, 2004, §7.1).
2.2 The boundary conditions of the flow and sediment
We formulate boundary conditions on the mean bed and on the free surface in terms of the turbulent mean velocity field , turbulent mean concentration , the fluid depth , and the turbulent mean pressure .
First formulate the boundary conditions for the turbulent flows.

On the mean bed, no fluid penetrating the ground requires ^{1}^{1}1This and the following boundary conditions are expressed in terms of ensemble mean quantities. Consequently, terms in the mean of the products of fluctuations might appear and a closure for them invoked (Cao, 2014, §2.2.6). We assume the closure that such products of fluctuations are negligible in the boundary conditions.
(8) 
Positing a slip law on the mean bed to account for a negligibly thin turbulent boundary layer gives
(9) (10) where the derivative and the constant matches open channel flow observations (Roberts et al., 2008, e.g.).

On the free surface (that is, on its turbulent mean position), the kinematic condition that no fluid crosses the free surface is
(11) 
Zero turbulent mean stress normal to the free surface indicates that on
(12) 
No turbulent mean, tangential stresses at the free surface indicates that on
(13) (14)
There are two boundary conditions for the suspended sediment.

On the free surface, the sediment flux normal to the surface is zero, which requires
(15) 
On the mean bed, the upward net flux across the mean bed comes from the entrainment due to the fluid turbulence, that is
(16) where, following the work of van Rijn (1984), the equilibrium reference concentration is approximated in terms of shear velocity and mean particle size as
(17)
2.3 Embed the physical problem in a family of problems
In order to provide theoretical support for the model redction, we embed the surface conditions (13)–(14) in a family of conditions that modify the tangential stresses to have an artificial forcing proportional to the square of the local, free surface, velocity:
(18)  
(19) 
When evaluated at parameter these artificial righthand side becomes zero and the artificial surface conditions (18)–(19) reduce to the physical surface condition (13)–(14).
However, when the artificial parameter , and when the mean slope and the lateral derivatives are negligible (), then the boundary conditions (18)–(19) reduce to and . In this case, two neutral modes of the dynamics are the lateral shear flow . ^{2}^{2}2The Euler parameter of a toy problem suggests introducing a factor into the lefthand side of the tangential stress boundary conditions (13)–(14) in order to improve convergence in the parameter when evaluated at the physically relevant . For the moment, this work omits such a factor.
Analogously we embed the sediment boundary conditions (15) and (16) in the family
(20)  
(21) 
Upon setting the embedded parameter , the boundary conditions (20)–(21) recover the original physical sediment boundary conditions (15)–(16). When the embedded parameter , boundary conditions (20)–(21) are part of an artificial problem which is used to find a slow manifold in the physical system. The extra term in equations (20)–(21) ensures conservation to errors . Without such term, the model only conserves sediment to errors . The reason for implementing such high order errors is that the information about different falling velocities comes into the model in terms.
When the artificial parameter , the lateral gravity and lateral derivatives are negligible (), and the falling velocity , boundary condition (20) requires the concentration on the free surface, and boundary condition (21) indicates
Together with the sediment pde (5), these imply a neutral mode of the sediment dynamics is when artificial parameter .
Conservation of fluid provides a third neutral mode in the dynamics. Thus, when , a four parameter subspace of equilibria exists corresponding to some uniform lateral shear, turbulent mean, flow, , some turbulent mean concentration , on a fluid of any constant fluid depth . For large enough lateral length scales, these equilibria occur independently at each location and (Roberts, 1988, 2008, e.g.) and hence the subspace of equilibria are in effect parameterised by , , and .
Provided we can treat lateral derivatives and as perturbing influences, that is provided solutions vary slowly enough in and , centre manifold theorems (Roberts, 1988; Chicone, 2006, e.g.) assure us of three vitally important properties:

this subspace of equilibria are perturbed to a slow manifold, whereon the evolution is slow, that exists for a finite range of gradients and , and parameters , , and , and which may be parameterised by the depthaveraged lateral velocities and , the depthaveraged concentration , and the local fluid depth ;

the slow manifold attracts solutions from all nearby initial conditions provided the spectrum of the linearised dynamics is suitable;

and that a formal power series in the parameters , , , and gradients and approximate the slow manifold to the same order of error as the order of the residuals of the governing differential equations.
That is, the theorems support the existence, emergence, and construction of slow manifold models such as (1). ^{3}^{3}3This general type of argument has recently been made rigorous in one lateral dimension by Roberts (2013). This support occurs in a finite domain in parameter space, and we assume the finite domain is big enough to include interesting cases of the physically relevant and finite slope and falling velocity .
2.4 Linear dynamics of the system
The linear dynamics of the system support the application of centre manifold theory. For the flat bed of , and with , the base problem (2)–(21) has the equilibrium of a shear flow which is, in terms of the stretched vertical coordinate ,
(22)  
(23) 
where and are the depthaveraged lateral velocities, is the depthaveraged concentration, and is the mean speed. Environmental turbulent flows have eddy viscosity variations. In this linear analysis, we assume the eddy viscosity is effectively constant.
Then we consider the dynamics of the pdes (2)–(5) linearised in the small perturbation fields about each of these equilibria. Because environmental turbulent fluids have very large lateral scales compared with the depth the lateral variations are very slow. As the lateral variations vary slowly they do not affect the dominant linear process. Thus we also treat the lateral derivatives and as ‘small operators’ in this linearisation. Thus, the pdes (2)–(5) and the boundary conditions (8)–(21), in the linearisation that effectively , result in the linear problem
(24a)  
(24b)  
(24c)  
(24d)  
(24e)  
(24f)  
(24g)  
(24h)  
(24i)  
(24j)  
(24k)  
(24l)  
(24m)  
(24n) 
Equations (24a) and (24f) determine there is no vertical velocity, . Equation (24i) implies the free surface perturbation , which corresponds to the freedom already in (22) so without loss of generality we here set . Equations (24d) and (24j) implies no change to the hydrstatic pressure, . pdes (24b), (24c) and (24e) indicate fields , and have solutions in the form
(25) 
where is a nondimensional vertical wavenumber. Substitute these solution forms (25) into pdes (24b), (24c) and (24e), and obtain . Substituting the solution forms (25) into boundary conditions (24k)–(24n) leads to two separate conditions for the velocity fields and concentration fields, respectively
(26) 
The characteristic equations (26) have the nonzero wavenumbers , which implies the leading nonzero eigenvalue . In addition to these negative eigenvalues, there are zero eigenvalues corresponding to the freedom to vary the fluid depth , depthaveraged velocities and and depthaveraged concentration . Thus, there is a spectral gap between the eigenvalues and . Centre manifold theory (Roberts, 1988; Potzsche & Rasmussen, 2006, e.g.) then supports the existence and emergence of a slow manifold of large lateral scale in the three dimensional turbulent fluid and sediment system.
2.5 Reduced model of the flow and sediment dynamics
This section focusses on interpreting the slow manifold of the leading order model of the turbulent flow and suspended sediment.
Instead of depthaveraging equations, the centre manifold approximation theorem underlies an iterative construction of a slow manifold that resolves the turbulent and sediment interactions within the fluid layer. The computer algebra program listed in Appendix A constructs the slow manifold of the turbulent flow and sediment system. The program derives evolution equations for the water depth , the depthaveraged lateral velocities and , and the depthaveraged concentration .
The order of error in the construction is phased in terms of the small parameters. Here the small parameters are the lateral derivatives and , the small mean slope , the falling velocity , and the artificial parameters and . Generally we report results to errors for some prescribed exponent (where, for example, a term with factor is if ). Further, the artificial small parameter is introduced to integrate the sediment dynamics into the theoretical support for a slow manifold. In order to ensure the sediment dynamics are accurate we construct the slow manifold to higher orders in the parameter .
Cao (2014) [§2.4.3] showed that coefficients of power series in converge quickly in the turbulent fluid flow systems. Computations with the sediment equations indicate that the dependence upon also converges reasonably quickly. Truncating to errors , the computer algebra in Appendix A derives the evolution of the depthaveraged concentration :
(27a)  
(27b)  
These power series converge quickly enough in to reasonably evaluate at to give accurate coefficients in the evolution equations.
Executing the computer algebra in Appendix A and evaluating at leads to the following evolution equations for the turbulent flow and sediment system in terms of the depth , depthaveraged lateral velocities and and depthaveraged concentration . The equations here are complicated due to the methodology systematically resolving all the intricate microscale physical interactions.
(28a)  
(28b)  
(28c)  
(28d) 
Equation (28a) is a direct consequence of conservation of fluid. The momentum equations (28b)–(28c) include the effects of drag , advection, such as , turbulent dissipation, gravitational forcing , and pressure gradients established by the suspended sediment . The sediment concentration equation (28d) contains the equilibration of vertical sediment distribution due to terms such as , advection such as , and turbulent dispersion effects.
Although equations (28a)–(28d) are expressed in terms of depthaveraged lateral velocities and depthaveraged concentration, they are derived not by depthaveraging, but instead by systematically accounting for interaction between vertical profiles of the velocity and concentration, the stress, bed drag, lateral space variations and bed topography. The form of coefficients in equations (28a)–(28d) are supported by dynamical systems theory: the detail in the equations reflects that a slow manifold is in principle composed of exact solutions of the Smagorinski dynamics (3) and convectiondiffusion equation (5), and hence accounts for all interactions up to a given order of analysis no matter how small the numerical coefficient in the interactions.
The sediment model (28d) includes all the dominated terms in established modelling (Wu, 2004; Duan, 2004; Duan & Nanda, 2006, e.g.). For example, Duan (2004) derived the following depthaveraged advectiondiffusion equation of suspended sediment:
(29) 
which consists of the effects of vertical distribution , advection and , and dispersion and . But the derived model (28d) also includes more subtle effects, which could be important for suspended sediment in complex flow regimes. The model (28d) further gives modifications in presence of the ratio due to different distributions of sediment in the vertical for the different levels of turbulent mixing characterised by the mean flow speed . The coefficients in equation (28d) are a little different to the established model (29). Take the advection term for example, the established model (29) has the coefficient , whereas our coefficient of such terms is . Physically, a higher falling velocity means sediment concentrates more near bed where the mean advection velocity is lower, and hence net transport will be slower. Our model (28d) has smaller dispersion coefficent compared with the model (29).
3 Crosssectional structures of the flow and sediment
The centre manifold approach does not impose a specific crosssectional velocity distribution as done by other methods, instead our approach empowers the sediment and Smagorinski turbulent equations to determine the crosssectional structures in outofequilibrium dynamics. Recall that the locally stretched vertical coordinate . This section concentrates on the vertical distribution of the lateral velocity and concentration in steady flow, and compare our prediction of the lateral velocity with analogous published experimental data (Schultz & Flack, 2007, 2013; Celik & Rodi, 1988, e.g.).
3.1 Distribution of the suspended sediment
We concentrate on the vertical distribution of the suspended sediment in the slow manifold of the system. Computer algebra in Appendix A derives the physical field of sediment concentration in terms of the depth , depthaveraged velocities and , depthaveraged concentration and stretched local normal coordinate on the slow manifold, evaluating at :
(30a)  
(30b)  
(30c)  
(30d)  
(30e)  
(30f)  
(30g)  
(30h)  
(30i)  
(30j)  
(30k)  
(30l)  
The vertical sediment distribution (30) describes the loworder shape of the slow manifold in state space. Physically this equation describes the details of the suspended sediment concentration in outofequilibrium flow. The terms in equation (30) have physical interpretations. For example, the line (30a) approximates the mean concentration, together with the lines (30b)–(30e), which forms the basic distribution of the concentration in the vertical in the presence of the depthaveraged concentration , the equilibrium bed concentration , the falling velocity , and the mean fluid speed . The lines (30f)–(30i) describe the effect by advection on the vertical distribution of the concentration. Lines (30j)–(30k) describe modifications due to the change of the bed topography. The line (30l) describes modifications due to lateral flow down a slope in its affect on the vertical distribution of the suspended sediment.
Now we explore the distribution of the suspended sediment in steady flows. Consider the turbulent flow of constant depth with suspended sediment flowing on a flat mean bed of constant slope ; that is, the mean bed . We consider the concentration fraction is small, . In this regime, in (28b)–(28c) the terms and are negligible. Then the evolution equations (28a)–(28c) predicts the equilibrium and , so the mean speed