Modelling suspended sediment in environmental turbulent fluids

Modelling suspended sediment in environmental turbulent fluids

Meng Cao School of Mathematical Sciences, University of Adelaide, South Australia 5005. or    A. J. Roberts School of Mathematical Sciences, University of Adelaide, South Australia 5005.

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, depth-averaged 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.

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 depth-averaging, resolves out-of-equilibrium 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 depth-averaging flow and sediment equations (Wu et al., 2000; Pittaluga & Seminara, 2003, e.g.). We explore the implications of changing the theoretical base from depth-averaging to a slow manifold of the turbulent Smagorinski large eddy closure.

Figure 1: Scene of Yellow River turning in Shilou, China’s Shanxi ( This river carries vast amounts of suspended sediment.

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 out-of-equilibrium dynamics of the coupled turbulent sediment system. We choose to parametrise the slow manifold in terms of emergent depth-averaged quantities. The evolution of these depth-averaged 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.12.2 detail equations of the Reynolds-averaged Navier–Stokes pdes, the advection-diffusion 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 non-dimensional suspended sediment model of the horizontal evolution of the depth , the lateral depth-averaged velocities  and  and depth-averaged concentration :


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, self-advection, 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 depth-averaged lateral velocities and concentration, they are derived not by depth-averaging, 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 out-of-equilibrium 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 advection-diffusion 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

Figure 2: This diagram depicts the suspended sediment in turbulent flow with  coordinate system. The fluid of depth  flows down the sloped bed at the turbulent mean velocity . The turbulent mean concentration of the suspended sediment is . Denote the mean bed , the free surface  , and the gravity is .

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 Reynolds-averaged continuity and momentum equations,


and for the suspended sediment is the advection-diffusion equation,


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


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


where the magnitude of the second invariant of the strain-rate 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.)


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 111This 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.

  • Positing a slip law on the mean bed to account for a negligibly thin turbulent boundary layer gives


    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

  • Zero turbulent mean stress normal to the free surface indicates that on

  • No turbulent mean, tangential stresses at the free surface indicates that on


There are two boundary conditions for the suspended sediment.

  • On the free surface, the sediment flux normal to the surface is zero, which requires

  • On the mean bed, the upward net flux across the mean bed comes from the entrainment due to the fluid turbulence, that is


    where, following the work of van Rijn (1984), the equilibrium reference concentration is approximated in terms of shear velocity  and mean particle size  as


The nondimensional pdes (2)–(5), together with boundary conditions (8)–(16) govern the dynamics of the turbulent flow and suspended sediment.

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:


When evaluated at parameter these artificial right-hand 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  . 222The Euler parameter of a toy problem suggests introducing a factor into the left-hand 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


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 depth-averaged lateral velocities  and , the depth-averaged 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). 333This 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 ,


where  and  are the depth-averaged lateral velocities, is the depth-averaged 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


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


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


The characteristic equations (26) have the non-zero wavenumbers , which implies the leading non-zero eigenvalue . In addition to these negative eigenvalues, there are zero eigenvalues corresponding to the freedom to vary the fluid depth , depth-averaged velocities  and  and depth-averaged 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 depth-averaging 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 depth-averaged lateral velocities  and , and the depth-averaged 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 depth-averaged concentration :


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 , depth-averaged lateral velocities and and depth-averaged concentration . The equations here are complicated due to the methodology systematically resolving all the intricate microscale physical interactions.


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 depth-averaged lateral velocities and depth-averaged concentration, they are derived not by depth-averaging, 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 convection-diffusion 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 depth-averaged advection-diffusion equation of suspended sediment:


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 Cross-sectional structures of the flow and sediment

The centre manifold approach does not impose a specific cross-sectional velocity distribution as done by other methods, instead our approach empowers the sediment and Smagorinski turbulent equations to determine the cross-sectional structures in out-of-equilibrium 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 , depth-averaged velocities  and , depth-averaged concentration  and stretched local normal coordinate on the slow manifold, evaluating at :


The vertical sediment distribution (30) describes the low-order shape of the slow manifold in state space. Physically this equation describes the details of the suspended sediment concentration in out-of-equilibrium 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 depth-averaged 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