Incorporating grain-scale processes in macroscopic sediment transport models

Sediment transport simulations face the challenge of accounting for vastly different scales in space and time that cannot be tackled by a unifying approach. Instead, processes are subdivided into a microscale at the particle level, a mesoscale of a large finite number of particles, and a macroscale that computes the sediment motion by means of advection–diffusion equations. The different processes occurring at different scales are simulated using different computational approaches. However, modeling sediment transport at multiple scales with high fidelity requires proper closure arguments that interconnect the different processes. Ultimately, we will need efficient macroscale models that can readily be utilized for engineering practices covering, e.g., entire river reaches or even estuaries. In recent years, highly resolved simulations have become a valuable tool to provide these closure arguments for sediment transport models on the continuum scale. In this paper, we will review the most relevant approaches to simulate sediment transport at different scales and discuss the perspectives of four most promising modeling techniques that can help to improve sediment transport modeling. On the grain scale, these enhancements include the impact of mechanical properties of cohesion and biocohesion as well as the shape of non-spherical sediment grains on fluid–particle and particle–particle interactions. On larger scales, we review constitutive equations for the macroscopic rheological behavior of sediment beds that may decouple the relevant scales for fluid and sediment motion. Furthermore, we discuss machine learning strategies as an efficient means to derive scaling arguments across multiple scales.


Introduction
Understanding sediment transport in rivers, lakes, and along the ocean floor is a key to sustainable management of open water bodies and aquatic ecosystems [1][2][3][4]. Owing to the increased computational power, numerical simulations of sediment transport have become a valid alternative to physical modeling in hydraulic laboratories or even at field scale [5]. This is especially true when it comes to the task of predicting the outcome of different scenarios after human intervention [6][7][8]. The most prominent examples comprise turbulent openchannel flows in fluvial hydraulics and the resulting morphodynamic evolution of river beds [9], turbidity currents as a particle-laden and gravity-driven flow down a slope [10,11], debris flows and mud slides [12,13], jets that lead to local scour of sediment beds at hydraulic structures [14,15], and coastal and estuarine morphodynamics that are influenced by low frequency oscillatory flows [4,16]. Owing to the different processes encountered at multiple scales in environmental fluid dynamics, the development of reliable predictive tools to simulate sediment transport phenomena across the various scales has emerged as a highly active and dynamic discipline in sediment transport research. Consequently, different processes have been classified depending on their scales in space and time [17], which allows to tackle these processes individually by using different computational approaches (cf. Fig. 1): (i) the microscale of individual particle-particle interactions can be B. Vowinckel (B) Leichtweiß-Institut für Wasserbau, TU Braunschweig, 38106 Braunschweig, Germany e-mail: b.vowinckel@tu-braunschweig.de Fig. 1 Concept of upscaling in sediment transport simulations after Seyed-Ahmadi and Wachs [18]: the continuum model shows a turbidity current propagating along the x-direction (from Zhao et al. [19]) represented by particle-resolved direct numerical simulations (DNS) that compute the motion of every grain by solving the Newton-Euler equations of motion and account for the fluid forces by resolving the flow field around it, (ii) the mesoscale that involves millions of particles is typically represented by particles smaller than the grid size of the fluid solver, so that grains are simplified as mass points and the relevant forces that govern their motion are computed by empirical correlations, and (iii) the macroscale can cover several kilometers, and the corresponding computational approach can no longer resolve the motion of individual grains. Instead, sediment transport is dealt with by continuum models of coupled transport equations. While these computational approaches can be sorted in terms of fidelity, they are also vastly different in terms of computational costs, with particle-resolved DNS being by far the most expensive.
Despite the high computational costs, particle resolved simulations are an important tool to identify the physical processes that remain relevant on larger scales and to provide a statistical description of these processes. These descriptions in turn can be used as constitutive equations for sediment transport models across scales. This means, however, that highly resolved simulations need to incorporate the most important physical effects, whereas limitations in computational capacity require idealizations for the description of the fluid-particle interaction. As two prime examples, many highly resolved studies have treated sediment as cohesionless, perfectly spherical grains [20,21], and it remains unclear, so far, as to how those assumptions distort the knowledge presently available at the microscale. As another trade-off, the more efficient point-particle and continuum models rely on empirical closures that inevitably are built on certain assumptions that may not hold for different situations. For example, while particle-resolved DNS resolves the fluid drag acting on a particle, larger-scale models typically account for fluid forces by using a drag coefficient that is derived from the fluid drag on a single sphere in an undisturbed flow, which is obviously not a valid assumption when dealing with sediment transport.
The present review addresses these issues. We will first review popular approaches to model sediment transport at different scales in Sect. 2 before we proceed with four active fields of research in Sect. 3 that bear a lot of potential for future enhancements of sediment transport modeling. It is clear that not all of the methods currently used can be discussed in these Sections. Therefore, the selection can only cover those examples that are considered to be the most promising, which inevitably may involve a slight bias from the previous work of the author [22]. On the microscale, more research on the role of cohesive inter-particle effects including the impact of biocohesion and the impact of non-spherical shape of sediment grains is needed to enhance the realism of particle-resolved simulations even further. On larger scales, the description of the rheological behavior of entire sediment beds needs suitable constitutive equations for meso-and macroscale sediment transport models. The rheological behavior can also lead to a decoupling of the relevant scales for the fluid and the morphology of the sediment bed, which may require additional considerations of the momentum balance for the fluid-particle mixture with regard to turbulence modeling. Finally, machine learning strategies can serve as efficient algorithms for upscaling relevant effects from the micro-to the meso-and macroscale. With the increasing computational power, all of these areas have recently seen a big push, since highly resolved numerical simulations can now provide the data needed to derive constitutive equations for larger-scale continuum models, but more work will be required in the future to exploit their full potential. It can be expected that the outcome from each of these fields will provide new fundamental insights into the governing mechanisms of sediment transport to transform our current understanding of sediment transport modeling. Consequently, this review is intended for communities using both highly resolved and macroscopic continuum-type sediment transport

Simulation strategies for sediment transport
Highly resolved numerical simulations of sediment transport deal with particles that are transported by an incompressible, viscous carrier fluid. Traditionally, this is done by solving the continuity equation, given by and the Navier-Stokes equations to compute the fluid motion: where u f denotes the fluid velocity, t is time, and ρ f indicates the fluid density. The fluid stress tensor is given by , where p f represents the fluid pressure, I is the identity matrix, and η f the dynamic fluid viscosity. The right-hand side also includes f bf as the sum of all body forces that may be acting on the fluid. These body forces can comprise a pressure gradient f p that drives a channel flow, a drag force f d that can potentially link particle motion to fluid motion, and a buoyancy force f b that can be caused by temperature, salinity or sediment concentration differences.
It is convenient to discretize equations Eqs. (1) and (2) on Cartesian grids. Depending on the spatial resolution of these grids, choices have to be made regarding the modeling of subgrid-scale turbulence. A proper DNS requires all relevant length scales of the flow to be resolved by the smallest grid cell size, and, hence, no subgrid-scale model is needed. On the other hand, a coarser grid may only resolve large eddies so that an excess diffusion, i.e., a turbulent eddy viscosity, enters the viscous term of Eq. (2) that needs to be modeled to account for the unresolved subgrid processes of the turbulent fluid motion. Depending on the spatial resolution and turbulence model, these types of simulations are either termed large eddy simulations (LES), where the largest, energy-containing eddies are still resolved by the fluid grid, or Reynolds-averaged Navier-Stokes (RANS) simulations or URANS for simulations of unsteady phenomena at larger scale, respectively. Since the focus of this review is the fluid-particle coupling and all three simulation techniques rely on the momentum balance expressed by the Navier-Stokes equations, we will not review the different turbulence models, but focus on the fluid-particle interaction. It is important to note, however, that the grid resolution will pose requirements for the proper representations of sediment-laden flows.
In recent years, much progress has been made regarding the modeling of flow-sediment interactions. There are different ways to run simulations of sediment transport phenomena depending on the geophysical phenomenon under consideration. The choice depends on three important dimensionless numbers, namely local sediment concentration φ as the part of a given control volume that is occupied by solids, the Reynolds number, and the Stokes number. The Reynolds number emerges from dimensional analysis by applying the Buckingham π-theorem to Eq. (2) [23]. Analogously, the Stokes number is the consequence of writing the equation of motion for a single-sediment grain as Newton's second law in a dimensionless form. The Stokes number St t p /t f as the ratio of the characteristic timescales of the particle (t p ) and the fluid motion (t f ), hence, is a measure to quantify the ability of a suspended particle to follow the fluid motion based on its inertia. Choosing the particle diameter and the maximum particle velocity as the characteristic length scale and velocity, the particle Reynolds number Re p max |u p − u f | D p /ν f , yields a description on the complexity of the flow field acting on the particle. Here, u p is the particle velocity vector, D p the particle diameter, and ν f η f /ρ f the kinematic viscosity. For wall-bounded flows, the flow becomes anisotropic in wall-normal direction. In this case, the smallest scales are characterized by the friction velocity u τ τ w /ρ f and the viscous length scale l τ ν f /u τ , where τ w is the wall shear stress [24], which reshapes the particle Reynolds number into a shear velocity Reynolds number D + D p u τ /ν f . Consequently, the particle behavior and fluid-particle coupling expressed by St and Re p or D + pose restrictions on the choice of the numerical model (cf. Fig. 2). This framework can also be recast in terms of a modified Shields diagram [25], where particle inertia is quantified by the Galilei number G D p g D p /ν f and the competition of destabilizing fluid shear stresses and the stabilizing weight of the particle is described by the Shields parameter θ u 2 τ /(g D p ) [26]. Here g (ρ s − 1)g is the specific gravity of the sediment, g is the gravitational acceleration, and ρ s ρ p /ρ f is the density of the solid relative to the fluid.
For St 1, Re p 1, and a small particle volume fraction φ ≤ O(1%), i.e., a dilute system with very small and light particles suspended in a flow, it suffices to solve Eqs. (1) , (2) on an Eulerian grid, and treat the sediment as a concentration field, whose spatiotemporal distribution is computed by an advection-diffusion transport equation on the same Eulerian grid. It is, hence, well suited to investigate larger-scale problems, such as dilute turbidity currents (Fig. 2) [2]. Based on the descriptions of the fluid and the particle phase, this method is commonly referred to as the Euler-Euler (EE) approach. According to Sommerfeld [27], the EE approach implies that a control volume is occupied by both phases, i.e., the fluid and the particle phase, at the same time, and their spatial and temporal evolution can be computed by applying certain interconnected conservation laws. In this spirit, we introduce the EE approach as an overarching framework that can be further subdivided into single-phase and two-phase flows. For the former, the particle phase is treated with an advection-diffusion transport equation of a continuous concentration field of a dilute particle phase that introduces a buoyancy term in the Navier-Stokes equations. For the latter, the particle motion is computed by a partial differential equation of the same form as the Navier-Stokes equation, so that a momentum exchange between the two phases can be described with better accuracy. Extending these methods to denser suspensions and particles of finite size is by no means trivial, but it shows potential for future uses given that appropriate closure laws can be found to describe the relevant physical phenomena of fluid-particle and particle-particle interaction involved in sediment transport.
For St ≥ 1, the trajectories of all the individual particles that are transported by the flow need to be computed in a Lagrangian sense [17,28]. This implies that the equation of motion is solved for every single particle. The fluid-particle coupling then mostly depends on the particle Reynolds number Re p . Particles smaller than the grid resolution may be treated as mass points, whereas the flow around larger particles needs to be geometrically resolved. Owing to the strategy of using a fixed Eulerian grid to numerically solve for the fluid motion and the Lagrangian approach to compute particle trajectories, this framework is commonly referred to as the Euler-Lagrange (EL) approach. Both of these approaches (EE and EL) will be discussed in detail in the following subsections. If, furthermore, φ ≥ 0.1%, particle-particle interactions due to collisions and contact can no longer be ignored [29], which calls for a coupling of the computational fluid dynamics with collision modeling from discrete element methods (CFD-DEM) as will be detailed in the following Sections.

Euler-Lagrange framework
The Euler-Lagrange approach solves for the fluid motion, and the sediment is described in a Lagrangian sense by computing the movement for every single grain by the Newton-Euler equations of motion for the translational velocity u p (u p , v p , w p ) T , and the angular velocity ω p (ω p,x , ω p,y , ω p,z ) T , Fig. 3 Different types of coupling schemes for particle-laden flows of different volume fractions. The limits that separate the different schemes were taken from [29] and [30] of spherical particles with a uniform density. Here, m p is the mass of particle p, F h, p the hydrodynamic force, F g, p the force due to gravity, I p is the moment of inertia of a spherical particle, T h, p is the torque due to hydrodynamic forces, and F c, p and T c, p are the forces and torques due to particle collisions, respectively. The computational approaches that make use of Eq. (2) and the particle equations of motion Eq. (3) and Eq. (4) can be categorized by their coupling [17] (Fig. 3): a one-way coupled system is obtained if particles are solely driven by fluid drag and gravity, which yields f d 0, F c, p 0, and T c, p 0; a two-way coupled system accounts for the feedback of particle drag on the fluid motion, which requires proper formulations for f d as will be detailed below; as a less common case, a three-way coupled system, where particles are driven by a flow and interact by particle collision and contact but do not modify the underlying flow field, is especially useful to investigate the formation of clusters and the agglomeration/flocculation of small suspended particles [30,31], but less ideal to investigate the hydro-morphological effects of sediment beds on the flow; and finally, a fully or four-way coupled system is realized by accounting for the momentum exchange between the fluid and the particle phase and computing particle-particle interactions of contact and collision. Note that there exists some ambiguity about the definition of three-way coupled systems, as they are sometimes defined as suspensions where particle flow disturbances affect the motion of nearby particles [32]. Since such a situation may be difficult to obtain without frequent particle collisions [27], this type of particle-particle interaction is defined here as another form of two-way coupled systems.
To realize particle-particle interactions for three-way and four-way coupling as defined in Fig. 3, discrete element methods (DEM) that are suitable for the given computational approach and vary by their degree of detail have been proposed. If particles are smaller than the smallest fluid scales that are resolved via the Eulerian grid to discretize the fluid domain, particles are considered as mass points, i.e., particles are described by their centers of mass and a virtual extent smaller than the grid size to attribute a particle mass [33]. On the other hand, particles that are larger than the smallest grid cell size of the Eulerian grid are considered to be spatially resolved. Both approaches make use of idealizing assumptions and bear different levels of sophistication as will be reviewed in the following Subsections.

Particle-resolved
For particle Reynolds numbers much larger than unity, the particle diameter exceeds the smallest spatial scale of the turbulent motion, i.e., the Kolmogorov scale, such that the flow around the particles has to be spatially resolved by the numerical grid. In recent years, a number of different methods have been developed to geometrically resolve mobile particles and to compute their motion in a viscous fluid flow. Prominent examples are in the fictitious domain the overset grid methods, PHYSALIS and the immersed boundary method (IBM), where the hydrodynamic forces and torques F h and T h are a direct result of the hydrodynamic stresses acting on the particle surface [34][35][36][37][38][39][40].
In the following, we will provide a description of DNS using the IBM with direct forcing, because this approach has become a popular method to simulate particle-laden flows [41][42][43][44]. The main advantage of the IBM is the use of Lagrangian marker points to discretize the particle surface. To this end, the particle surface is treated as a shell that is partitioned in a number of equal-sized cells as indicated by the orange marker points in Fig. 4. Efficient interpolation between the Lagrangian markers and the underlying Eulerian grid allows to compute the hydrodynamic forces and torques acting on the particle as the surface integral of the hydrodynamic stresses [41]: and where Γ p the fluid-particle interface, n is the outward-pointing normal on the interface Γ p , and r X l − x p is the position vector of the surface point X l with respect to the center of mass x p of the particle. Vice versa, the fluid at an arbitrary location x experiences drag due to particle motion via f d f ibm : where N tot is the number of particles in the domain, N l, p is the number of Lagrangian markers to discretize the surface of particle p, δ h is a Dirac delta function that spreads the forcing from the point on the particle surface to the Eulerian grid [45], U d l, p is the desired velocity at the particle surface, U l, p is the interpolated velocity computed prior to the forcing correction, and V l, p is a volume element of the thin shell surrounding the particle surface as represented by the Lagrangian marker points. Hence, f ibm acts as a correction to enforce the no-slip condition on the particle surface. It is also important to note that by using r X l, p − x p and I p const, the set of equations Eqs. (3)-Eq. (7) incorporate the assumption of rigid, spherical particles. The rigorous fluid-particle coupling by Lagrangian marker points allows for the efficient computation of particle motion by means of ordinary differential equations without remeshing the Eulerian grid for the fluid solver. In addition, dealing with denser particle systems requires suitable treatment of particle-particle interactions to obtain fully coupled flow conditions without loosing the advantages of the IBM. This is especially true for the situation of sediment transport, where the typical situation comprises the simultaneous collisions of fast moving particles and enduring contact of the grains within the sediment bed. Consequently, several sophisticated collision models have been proposed to account for particle-particle interactions addressing the particularities of the IBM coupling [46][47][48][49]. All of these studies have in common that the contact forces can be described by a set of spring-dashpot systems. This set comprises a repulsive and a frictional force acting in the normal and tangential direction, respectively. In addition to that, the tangential component is limited by the Coulomb friction law. The modeling framework is sketched in Fig. 5. While there are strong similarities in the contact modeling described in the references cited above, we present a brief summary of the approach proposed by Biegert et al. [49] to illustrate the general concept. The idea is to subdivide a particle collision into three phases: approach, contact, and rebound. During the approach phase, the gap size ζ n of the colliding objects becomes smaller than the grid cell size h. Hence, the work performed to squeeze the fluid out of this gap can no longer be captured by the IBM but has to be modeled by a lubrication force [50]: where the subscript pq indicates interactions of particle p with particle q, R eff R p R q /(R p + R q ) is the effective radius of the two interacting particles p and q, ζ min is a limiter that can be interpreted as the particles' surface roughness, and g n is the normal component of the relative velocity of the two colliding particles. During contact, a soft-sphere model represented by a nonlinear spring-dashpot model for the normal direction prevents particles from non-physical overlap, where k n and d n represent stiffness and damping coefficients according to the Hertzian contact theory [51] and g n is the normal component of the relative velocity. The coefficients k n and d n can be adaptively calibrated for every collision or contact to prescribe a desired restitution coefficient e −u out /u in [46], where u in and u out are the impact velocities of two particles or a particle and a wall right before and after the particle collision, respectively. The forces in the tangential direction are modeled by a linear spring-dashpot model capped by the Coulomb friction law as where k t and d t are stiffness and damping for the tangential force, g t is the relative tangential velocity, t is the unit vector pointing in tangential direction, μ f represents the coefficient of friction between the two surfaces, and ζ t is the tangential displacement integrated over the time interval for which the two particles are in contact [52], which has first been proposed for DEM of dry granular media [52,53]. Finally, as particles rebound, the lubrication model Eq. (8) is needed as long as ζ n < 2h until the full hydrodynamic force can be computed via the IBM forcing Eq. (7), so that the total force due to contact and collision entering Eq. (3) becomes: It is important to note that using a sub-stepping routine turned out to be a key ingredient for the fluid-particle coupling via the IBM [48,49]. By using sub-stepping, additional integration steps are introduced that are smaller than the original integration time step of the fluid solver. During these substeps, only equations Eq. (3) and Eq. (4) are solved by recomputing the collision forces, whereas the fluid forces acting on the particle are held constant during the duration of a time step of a fluid solver. This measure is needed to solve the inelastic collision by means of a soft-sphere model and to properly resolve the lubrication forces during approach and rebound in time.
A number of studies have recently been published that successfully applied the concept of fully coupled CFD-DEM outlined above to investigate sediment transport phenomena by means of particle-resolved DNS. Using the collision model of   [46], Vowinckel et al. [20,54,55] carried out simulations of bed-load transport of monodisperse spherical particles in a turbulent open-channel flow with a shear velocity Reynolds number of D + u τ D p /ν f in the transitionally rough regime (D + 21) [56]. In these studies, different modes of transport and sediment pattern evolution were investigated. The dataset presented in [20] was statistically analyzed for incipient motion events of individual particles and gave a clear picture of the coherent structures that are responsible for particle dislodgment [57].
Additionally, work has been devoted to the analysis of particle transport in laminar shear flows [58] as well as the trajectories and turbulence modulation of small particles with D + 5.6 in the hydraulically smooth regime [59]. It is clear, however, that a large number of particles will give rise to the development and evolution of ripples, which will increase the hydraulic roughness beyond the hydraulically smooth regime. This has been the scope of the studies by [21,60], who performed phase-resolved simulations of a unidirectional, steady turbulent open channel flow laden with more than 1 million spherical, inertial particles with D + 12, and [61,62], who carried out similar simulations under oscillatory flows. The computational domains of these studies were large enough to produce a sequence of ripples under both flow conditions.

Point particles
For particles smaller than the Eulerian grid cells, the particle surface is no longer geometrically resolved, but it suffices to treat the particles as mass points with some virtual extent that remains smaller than the grid cell size [17]. Note that this requirement does not necessarily prescribe Re p < 1, since computational approaches employing subgrid-scale turbulence modeling such as LES can handle particles larger than the smallest fluid scales as long as they are still smaller than the grid cell size of the Eulerian grid [28]. Nevertheless, the fact that particles are not spatially resolved requires empirical correlations to compute the drag on the particles. These correlations are typically derived from a sphere settling under its own weight in a quiescent, viscous fluid [63][64][65]. The relevant effects comprising the hydrodynamic forces on the particles are the drag (F d, p ) and lift forces (F l, p ), effects of the pressure gradient from the undisturbed flow (F f, p ), the added mass (F a, p ),and the Basset history force (F b, p ) as summarized by [66]. It is obvious that this approach introduces several empirical closures to model the overall drag when the flow around the particle is no longer resolved. In particular, fluid drag force is usually approximated as being pressure-dominated: where A p is the projected area of a particle settling in a fluid and C d is a drag coefficient. However, it is well known that C d is Reynolds number dependent [67]. While there are well established relationships for individual particles, closures for the effective drag in a dense system are more difficult to derive [68,69]. Computing fluid drag via Eq. (3) in combination with Eq. (12) is commonly referred to as the Basset-Boussinesq-Oseen (BBO) equation. This equation and variants thereof have gained a wide popularity for Euler-Lagrange-type simulations of sediment transport in turbulent flows. Even though it was pointed out by Balachandar [28] that point particles are particularly useful for LES, since the coarser grid allows to accommodate particles larger than the smallest fluid scales, this choice bears challenges on its own. Recent studies of sediment transport compared results from DNS and LES in the same flow regime to find that the two simulation strategies do not provide identical results [70]. Selected examples of point-particle EL simulation campaigns that investigated sediment transport in turbulent flows are given in Table 1 and many more could be mentioned here. It should be noted that the BBO equation is not strictly bound to using those four components of hydrodynamic force (drag, lift, added mass, and Basset history force), but can be simplified or extended as indicated by the (...)-term in Eq. (12) depending on the physical system under consideration such as the consideration of the Magnus force in shear flows [71,72].
As a pioneering work, a relatively simple one-way coupled approach of the BBO equation that merely accounts for drag and lift was used in seminal numerical studies of particle-laden flows in homogeneous turbulent shear and open-channel flows [73,74]. Despite this simple one-way coupled modeling strategy,  [77] X X X X X Barati et al. [72] X X X X X X Zhang et al. [78] X X X Zhao et al. [31] X conclusions could be drawn on local particle concentrations and dispersion. On the other hand, two-way coupling is needed to investigate turbulence modulation by inertial particles [79][80][81]. For this approach, a body force for the hydrodynamic drag f d is introduced in Eq. (2) as the sum of the drag force on all particles within a integration cell of the computational grid. Studies focusing on sediment transport, however, typically involve heavy particles that are transported along the channel floor or even form a granular sediment bed. To this end, a four-way coupled scheme is needed to take particle contact and particle collisions into account. Two types of collision modeling are typically used: (i) hard-sphere models and (ii) soft-sphere models. The former method has been developed for point-particle simulations in the dilute regime [82]. The idea is to have instantaneous collisions and prescribe the particle velocity by the restitution coefficient. Using this method, simulations of saltating particles along a sediment bed were carried out [71,72,75]. This method is well suited for suspended particles colliding with rather large impact Stokes numbers, where u in is the impact velocity right before the collision. However, owing to the wide range of impact Stokes numbers for particle collisions and sustained contacts during sediment transport, particle interactions of entire mobile sediment beds are typically simulated by a soft-sphere model. This type of contact model is similar to the framework laid out in Sect. 2.1.1. It is based on a combination of spring-dashpot systems that reflect repulsion and friction in the normal and tangential direction, respectively. To this end, the model allows for a small overlap of the colliding particles, and the contact duration is slightly stretched in time to numerically resolve the particle interaction [76][77][78]. These types of simulations were able to reproduce bed-form evolution of ripples and dunes in unidirectional and oscillatory flows that were very similar to the results reported for particle-resolved DNS [21,62,83].

Single-phase flow model
Many environmental or geophysical flows involve suspensions with low volume concentrations and, hence, small density differences. If the particle size in terms of Re p and St is much smaller than unity, the particles respond rapidly to the fluid flow, and if, in addition, the particle concentration is small enough (O(1%)), the density difference of the suspension and the effect of particle-particle interaction becomes negligible for the fluid inertia [2,84,85], i.e., ρ f remains constant in Eq. (2). In this case, it can be assumed that particles are advected with the fluid and experience a settling velocity based on the Stokes drag force. The fluid-particle interaction entering Eq. (2) then becomes the buoyancy force: where u s is the Stokes settling velocity, Fig. 6 Concentration iso-surface of a turbidity current propagating down a slope [23] e g is the unit vector pointing into the direction of gravity, and c s is the concentration of the suspension. On the other hand, the spatiotemporal evolution of the sediment suspension can be described by an advection-diffusion transport equation, where α is the diffusion coefficient of the sediment and R represents sources and sinks due to resuspension and deposition, respectively. Since density differences are small, they are neglected in the Navier-Stokes equation (2), but enter the momentum balance of the fluid only via the buoyancy force Eq. (15). This strategy is commonly known as the Boussinesq approximation [86]. A detailed derivation including the non-dimensional set of equations can be found in, e.g., [23]. It has been shown that this simplified model is an useful approach to investigate turbidity currents such as underwater avalanches [2,19,[87][88][89] (Fig. 6) or the transport of fine sediment in oscillatory flows [90,91].
It is important to note that the distinction between using an Euler-Euler and Euler-Lagrange coupling is not strictly bound to the thresholds indicated in Fig. 2, but the parameter ranges can overlap. For example, a recent EL simulation by Shao et al. [92] produced very similar results compared to the EE simulations of Burns and Meiburg [93], which both confirmed results from linear-stability analysis [94].
While simulations of turbidity currents typically incorporate the assumption of a smooth wall that adsorbs sediment as it settles to the ground but does not allow for resuspension, this situation is not true for sediment transport, where erosion processes can add sediment mass to the concentration field as a boundary condition. Erosion and resuspension can also lead to the development of bed forms. This has been addressed in recent years, by extending the model with the Exner equation as the lower boundary condition for open-channel flows [95][96][97][98]: where h is the height of the sediment bed, is a diffusion term, and q is the volumetric sediment flux that reflects erosion, resuspension, and deposition. Another important assumption of the single-phase flow model presented by Eqs. (15)- (18) is that particles are assumed to be cohesionless and settle at the rate of an isolated spherical particle in an otherwise undisturbed environment. These considerations can be extended to cohesive sediments by applying formulae for the suspension and the deposited sediment to account for hindered settling [99] and consolidation [100], respectively. While there has been a number of studies that have taken these processes into account, the closures commonly used remain empirical [101,102].

Two-phase flow models
For obvious reasons, the assumption of a dilute suspension to solve Eq. (2) for the fluid motion is not always applicable for sediment transport modeling. A more comprehensive extension of the EE framework for larger St and the erosion/resuspension and deposition in denser granular systems has been derived in terms of twophase flow models [103][104][105]. These types of models are based on the Reynolds decomposition ξ ξ − ξ for a velocity component ξ of either the fluid or solid phase, whereξ is the ensemble average and ξ is its fluctuation obtained from the Favre averaging operator [106,107], to distinguish the fluid and the solid phase [108]. Here, the subscripts f and s indicate the fluid and solid phase, and ψ and φ 1 − ψ represent the respective local porosity and solid volume fraction averaged over the same time interval. The averaging operator Eq. (19) yields modified conservation laws for mass and momentum for the fluid: and for the solid phase: whereũ s is the velocity vector of the solid phase and p s is the granular pressure. Just like the single-phase momentum balance Eq. (2), the two-phase model balances the material derivative comprising local acceleration and convective transport with the internal forces originating from the pressure gradient and the shear stress as well as external body forces. Note, however, that these external forces now include the momentum exchange between the two phases I, which is a function of the local velocities and volume fraction, i.e.,ũ f ,ũ s , and φ, respectively, as well as the drag coefficient C d of the solid particles making up the solid phase. In addition, the shear stress tensors T f and T s can be modeled by means of an effective viscosity that accounts for effects from turbulent mixing on the subgrid-scale and particle-particle interactions. Since the Favre average Eq. (19) introduces filtered equations, the fluid stress can be modeled with closures from classical turbulence modeling such as RANS [109] or LES [110], whereas the particle stress component requires additional treatment regarding the rheological properties of the granular system as will be discussed further in Sect. 3.3 below. Furthermore, a few two-phase fluid models exist that are able to unify sedimentation and hindered settling [111][112][113]. It was shown in these references that the two-phase flow model allows for a more rigorous coupling of the fluid and the particle phases than the single-phase flow model. However, these examples still deal with a simplified one-dimensional sedimentation problem, i.e., cohesive sediment settling in a water column, and they inherit the empirical closures that were already pointed out as a shortcoming for the single-phase flow models in Sect. 2.2.1.

Future applications and enhancements
The computational frameworks presented in Sect. 2 are sorted from the treatment of small to large systems. In the following, we present four most promising examples to extend these frameworks towards a higher degree of realism. This includes the exploitation of EL-type simulations that can be enhanced to account for cohesive forces and the effects from non-spherical particles on the grain scale, machine learning techniques as an efficient means to transfer the knowledge gained on the microscale to larger-scale systems, and finally frameworks of constitutive equations that allow to close the systems of equations of EE-type simulations summarized in Sect. 2.2 above.

Cohesive and adhesive sediments
The particle-particle coupling by the contact modeling framework presented in Section 2.1 only holds for cohesionless sediment such as sand and gravel. However, many applications in environmental engineering involve smaller-sized grains, such as silt (2μm ≤ D p ≤ 63μm) and clay (D p < 2μm). Those types of sediments are commonly found in aquatic environments including larger rivers, estuaries, coastlines, and the ocean floor where the local flow rate decreases, which yields a decreasing transport capacity within the water body. For silt and clay, cohesive van der Waals (vdW) forces can dominate over gravitational forces and trigger particle aggregation or flocculation.
Larger-scale EE-type models that are used for the management of entire estuaries and harbors have to incorporate these processes for the transport of mud, which is a mixture of silt, clay, and organic matter. One such model was first proposed by Waeles et al. [6] for the Seine estuary by a number of advection/diffusion equations for different grain sizes. The model is designed to cope with domain sizes that have scales of several kilometers [7,114]. Since such a macroscopic modeling approach relies on a number of empirical coefficients that need to be recalibrated for every new physical scenario, more work will be needed to describe the implications of small-scale processes such as particle flocculation and settling as well as resuspension, so that suitable subgrid-scale models can be found for EE-type simulation strategies on the macroscale.
In this regard, EL simulations can provide the data needed, and, indeed, recent studies have indicated that this is possible. For example, point particle approaches have utilized the "hit and stick"-condition by evaluating the van der Waals forces, and variants thereof [115][116][117][118], where A H is the Hamaker constant, R eff the effective radius of two interaction particles of different size, and ζ 0 is the minimal separation distance. Point particle simulations can reproduce agglomeration phenomena, but they may suffer from the shortcomings of the empirical correlations that are used to compute the hydrodynamic drag on individual particles. This can be avoided by fully resolved EL simulations, albeit the increasing computational cost or smaller spatial scales, respectively [119][120][121].
Particle-resolved simulations also allow for contact modeling that leans more into the classical DLVO theory of Derjaguin and Landau [122] and Verwey and Overbeek [123] for colloidal particles. To this end, a cohesive range λ can be defined as a finite-sized shell surrounding the particle. This allows to define a cohesive force F coh as a function of the normal distance from the particle surface. If another object is within that cohesive range, attractive forces will pull the particle towards the other object. For example, Pednekar et al. [124] applied the DLVO theory by superposing the repulsive double-layer force with the attractive van der Waals force to obtain a total force that induces cohesiveness: where A r is a repulsive force scale and κ is the Debye length (Fig. 7a). A simplified version that cuts off the exponential tail of decaying cohesive forces with larger distances was proposed by Vowinckel et al. [120,121] to address the needs of particle resolved simulations: Note that this model is parameterized by a single parameter, the cohesive number Co max(||F coh ||) m p g (27) as the ratio of the maximum cohesive force at λ/2 and the weight of the particle. A sketch of this model is given in Fig. 7b and c. The highly resolved data from EL simulations allows to quantify the hindered settling of complex aggregates [99,120] as a means to parameterize the settling velocity in Eq. (17). The simulation results of polydisperse sediment settling in a quiescent fluid by Vowinckel et al. [120] confirmed the enhanced settling speeds of  [120]. The dashed red line in (a) demarcates the range λ that is described by the parabolic model by cutting off the exponentially decaying tail. (c) schematic of the cohesive range λ over which cohesive force acts during particle-particle interaction flocculated cohesive sediments. The highly resolved data of the grain trajectories presented in this study revealed that the enhanced settling speeds of cohesive sediments are caused by smaller grains trailing larger ones (Fig. 8).
Additionally, these types of simulations have the potential to quantify the effective normal stress [125] and the effective viscosity of the fluid-mud mixture [126] as cohesive sediment is transported, deposited, and consolidated over time, which are the closures needed for the models described in Sect. 2.2. The spatiotemporal evolution of cohesive aggregates can also help to improve existing and derive new population balance models to predict floc sizes and settling velocities, which has been a popular approach to parameterize advection/diffusion equations, that can reflect floc size distributions by employing monodisperse, two or three classes of sediment types [127][128][129]. A recent study by Zhao et al. [31] demonstrated that a three-way coupled system, i.e., neglecting the particle drag on the fluid flow field and modeling particle-particle interaction, is a very efficient means to investigate flocculation processes of cohesive particles. In this study, particles were driven in a steady, two-dimensional cellular flow field that was prescribed by the initial conditions for Taylor-Green vortices. This setup provided an efficient model to gather flocculation data for the parameterization of population balance equations that describe the growth of aggregates from initially dispersed grains of cohesive sediment.
It is important to note that the particle-resolved DNS discussed above [120,121] were meant to represent the behavior of macroscopic grains that correspond to an average diameter of 20μm. At this size, particles are small enough for cohesive forces to become important, i.e., smaller than 63μm, but large enough to behave as inertial particles that may be treated using the scheme described in Sect. 2.1.1. For smaller particles sizes, the particle Reynolds number becomes so small that the method of choice should be the point-particle approach. Indeed, it was shown by [31] for a grain size of 5μm that this is a reasonable approach even for one-way coupled simulations, where the impact of the particle motion on the fluid flow is neglected. For even smaller particle sizes, that is grains smaller than 2μm, one would have to include an additional randomly fluctuating force F b in Eq. (3) that introduces Brownian motion [130,131]. This force has the properties that the ensemble average is equal to zero at a given instant, but yields a net drift over time [132]. This type of modeling has successfully been applied in Stokesian dynamics simulations of spherical colloidal particles [124,133]. For these types of simulations, the hydrodynamic forces can be approximated by the Stokes drag for spherical particles F h −3π D p η f (u p − u f ). On the other hand, it is important to note that the assumption of a spherical particle does no longer apply to naturally occurring cohesive sediment with a grain diameter of less than 2μm. The granular material of this size consists of clay particles that are shaped as platelets. Hence, special treatment of non-spherical particles becomes necessary for this type of sediment.
While first results of EL simulations show the potential to investigate cohesive sediment dynamics, biofilms as another key component have often times been neglected, but may also induce agglomeration. Biofilms can grow on any substrate, but especially on fine grained sediments, because they are rich in nutrients that provide an excellent fertile environment for the growth of biofilms [134,135]. These biofilms can form a matrix of extracellular polymeric substances (EPS), which results in a glue-like coating of the sediment grains. This coating is highly adhesive and results in so-called biocohesion or biostabilization of the sediment, which can have a strong impact on the fluid-sediment interaction [136,137]. Even though the cohesion model of Vowinckel et al. [120] was derived for electrostatic interaction, an extension towards biocohesion seems to be feasible [121], but more work will be needed to quantify cohesive forces from actual experimental data and to account for the glue-like particle bonding.

Non-spherical particles
The computational frameworks described in Sect. 2 all have in common that they rely on the assumptions of spherical particles. Spherical particles have the advantage that they are uniform in shape and perfectly symmetric. This allows for the description of hydrodynamic drag on small particles via the analytical solution of Stokes drag in creeping flows, i.e., equation Eq. (16). At higher particle Reynolds numbers Re p , it was shown that empirical laws for fluid drag become applicable [33,67] that enter the BBO equation in terms of drag, lift, Basset and history force. In addition, the symmetric body makes it conceptually simpler to compute particle inertia and the particle-particle interaction. However, it is obvious that particles in natural environments are usually not perfect spheres but can have a variety of arbitrary shapes [138,139], so that the assumptions for spherical particles made in Sect. 2 above are no longer valid [140]. Hence, the impact of shape on the transport behavior will be a crucial issue in future studies to address the complexity of natural sediments.
One of the most prominent examples of experimental evidence that acknowledges this more complicated situation has been presented by Schmeeckle et al. [141], who investigated the collision behavior of natural sediments and compared it to spheres. Most notably, this study showed that the more complex shape of natural sediments will lead to lower restitution coefficients and, hence, lower rebound heights of non-spherical particles colliding with a high-impact velocity, i.e., the impact Stokes number St > 39 as defined by Eq. (14). Subsequently, four-way coupled LES-DEM have been performed with point particles and a simplified version of the BBO equation, where added mass, Basset and lift force are neglected [142]. The drag relation is corrected for non-spherical shapes by using a shape factor [143,144]. While these simplifications allowed for first simulations of sand particles, the idea of the shape factor has recently been challenged by results from highly resolved simulations [145,146]. Thus far, it remains an open question whether the description of the particle shape by one single parameter, e.g., sphericity or shape factor, is the way forward to account for this complex situation. Hence, more detailed investigation on the particle scale is needed to upscale the effect of more complex shapes on the larger-scale flow behavior and the different modes of bed load transport (saltating, sliding or rolling [147]).
One way to better approximate objects of more complex shape is the bonded-sphere approach, which has been applied to EL simulations of point particles [148,149]. In these studies, particles consist of several spheres that are rigidly bonded in one larger particle and the total size of that particle needs to exceed the grid cell size to resolve the resulting shear forces and torques. It is ultimately desirable, however, to use particle-resolved EL frameworks that directly compute fluid drag and particle-particle interactions acting on the natural sediment. Four-way coupled simulations of particle-resolved EL simulations, however, can no longer use the simplifying geometry of spherical bodies. On the one hand, the geometry requires special treatment to compute the particle motion. The complex shape of the particle needs to be discretized by Lagrangian markers, and the particle as well as its moment of inertia I p that appears in Eq. (4) needs a proper description. Equation (4) readily makes use of the assumption that particles are spherical and fully symmetric so that I p is isotropic and, hence, becomes a scalar, which allows for the simplification of Eq. (4). This simplification is no longer valid for ellipsoids, and the full tensor needs to be rotated according to the particle orientation. Here, quaternions are the method of choice to keep track of the particle orientation over time [140]. On the other hand, challenges arise for the computation of particle-particle interactions. Simulating a multitude of particles requires efficient contact detection [150] and the computation of local collision forces at contact points with non-uniform curvature [151]. In addition, the torque due to normal contact forces needs to be accounted for since the impulse is no longer directed through the center of mass.
One of the first studies to investigate a system of many non-spherical particles was carried out by Shardt and Derksen [152], who simulated the settling behavior of blood cells. The surface of blood cells was discretized by a random placement of markers with a characteristic distance that reflects the grid cell size of the fluid grid. The particle-particle interaction was dealt with in a minimalistic way. Contact points were determined by repeatedly checking the distance of surface markers that were organized in linked lists and contact forces merely comprise a repulsive force, whereas tangential and lubrication forces were neglected. In this study, the blood cells were found to align along the direction of gravity to minimize their drag forces during the settling process, which confirmed experimental observations.
A straightforward extension of a spherical particle is the shape of ellipsoidal objects. Ellipsoidal particles can be further subdivided depending on the ratio of their half axes a, b and c: (a < b and b  c), and (iv) ellipsoid (a b c). Options (ii) and (iii) are the most obvious extensions of the spherical geometry, because they are created by rotating an ellipse about one of its half axis. These types of particles are referred to as spheroids (cf. Fig. 9). An early study by Zastawny et al. [153] investigated the hydrodynamic drag and lift forces as well as torques for spheroids depending on the orientation and rotation to the flow direction. In this scenario, the particle was exposed to a uniform flow, and its center of mass was fixed in space and time. This study allowed to extend the BBO equation for neutrally buoyant spheroidal particles by providing valuable information for low concentrated suspensions where particle-particle interactions can be neglected. However, it is not straightforward to extend such a concept to denser systems that are contact-dominated, e.g., bed load transport or the rheology of dense suspensions. Moreover, sediment transport of gravel-sized particles (D p > 2 mm) often involves particle Reynolds numbers that are substantially larger than the viscous length scale, so that the flow around the particles has to be fully resolved [17] (cf. Sect. 2.1.1).
Consequently, a number of algorithms that are all based on the IBM have been developed to carry out particle-resolved DNS with spheroids [154,155] and the even more general shape of ellipsoids [156]. These shapes have the advantage that there will be only one contact point upon collision, which provides a straightforward extension of the DEMs that were originally proposed for spherical particles. Thus far, DNS that investigate the interaction of fully resolved non-spherical particles have mostly been dealing with suspensions of neutrally buoyant non-cohesive and non-Brownian spheroids [157][158][159] to report drag reduction by preferential alignment of the spheroidal particles (oblates and prolates) with the flow direction. These results confirmed the observations previously reported for point-particle approaches [140].
However, as already discussed in Sect. 2.1 proper DNS of sediment transport phenomena need to account for the more challenging situation of simultaneously occurring high-impact collision and sustained contact in dense granular systems. Recently, promising studies have appeared that addressed incipient motion in laminar and turbulent flows by means of particle-resolved DNS [160] and the effect of particle shape on sediment transport [146]. Both studies reported that while the particle shape does not seem to have an impact on the mean streamwise velocity, resuspended spherical particles jump higher and are transported with higher speeds thereby confirming the experimental observations of Schmeeckle et al. [141]. When in contact with the sediment bed, spherical particles are more likely to roll across the rough sediment bed, while the non-spherical particles tend to slide along the flow direction. On the other hand, the lower traveling speeds of ellipsoidal particles impose a larger hydraulic resistance, so that the turbulent fluctuations, i.e., the Reynolds stresses, are enhanced (Fig. 10). These results suggest that, indeed, the particle shape has a strong impact on the mode of sediment transport, such as rolling versus sliding [147], and the resulting flow behavior.
These considerations of particle-particle interactions can be further extended to particles of arbitrary shape, where multiple contact points per collision can occur, which would further complicate the picture [161]. Similar considerations should hold for the fluid drag on more irregular particle shapes as described, e.g., by [145] or the considerations of clay platelets discussed in Sect. 3.1 above. It will be of great importance in future upscaling strategies to reflect these properties in scaling laws and constitutive equations that are still to be derived for point-particle and EE-type simulations.

Rheology modeling
The modeling framework comprising the conservation laws for mass Eq. (1) and momentum Eq. (2) in conjunction with the scalar transport and Boussinesq approximation Eqs. (15)-(17) is a well-posed problem, and has been used with great success to model turbidity currents [2,84,85,162]. However, this framework is limited to the assumption of a dilute suspension with only small density difference compared to the clear water ambient fluid. An extension to denser systems as laid out by the two-phase fluid model Eqs. (20)-(23), however, requires constitutive equations for the fluid-particle interaction to close the system of equations, which is an active field of research, because the equations need to reflect fluid and particle properties. Hence, we will focus on closure strategies for the sediment part of the two-phase flow model presented in Sect. 2.2.2 in the present Section. In particular, strategies are needed to describe the effective viscosity of the mixture, the granular pressure and shear stress tensor, and the exchange of momentum between the fluid and the particle phase.
For suspensions of inertial particles, it is well known that the resistance to deformation of the suspensions is no longer governed by the internal molecular viscous stresses alone, but also by the inter-particle forces due to contact and friction. This gives rise to the idea of an effective viscosity η eff that is larger than the clear-fluid viscosity but recovers the original viscosity of the liquid for the clear fluid, i.e., η eff → η f as φ → 0. At the dilute limit, it is possible to develop analytical solutions by deriving approximations of first and second order [163,164]. Unfortunately, these analytical solutions are not applicable for dense suspensions that are typically encountered in sediment transport. Owing to the nature of sediment transport, the effective viscosity can be further subdivided in a shear and normal viscosity component defined as η s τ/(η fγ ) and η n p s /(η fγ ), respectively, where τ is the total shear that comprises fluid and particle stresses, andγ is the shear rate.
Empirical correlations for η s (φ) and η n (φ) have been proposed for experiments with volume-imposed rheometry to provide scaling with respect to the local volume fraction covering the dilute and the dense regime [165,166]. These correlations have two key features. First, they recover the clear-fluid viscosity for unladen flows, and, second, they prescribe η eff → ∞ as φ → φ c , where φ c represents the maximum critical packing fraction. At this limit, the particle dynamics enter a quasi-static regime, which is the so called jamming transition [167]. The effective viscosities η s and η n directly enter T s in Eq. (23). The normal viscosity also provides a means to prescribe the particle pressure p s as the total submerged weight within a given control volume. Unfortunately, these correlations involve empirically determined parameters. For example, Morris and Boulay [166] found φ c 0.68 for their experiments on the shear-induced migration of monodisperse particles; a value exceeding the randomly close packing. Consequently, it is clear that the empirical coefficients may have to be readjusted for different flow scenarios and, hence, lack general predictive capacity.
A step forward in this regard has been made by Boyer et al. [168] by transferring considerations of dry granular material [169] to viscous suspensions in a pressure-imposed setup, where the particle packing is free to dilate under shear. To this end, the dimensionless inertial number I D pγ ρ p / p s that describes the competition of strain and rearrangement of a sheared dry granular material was replaced by the viscous number J η fγ /p s to formulate empirical scaling laws of the macroscopic friction factor μ(J ) and the local volume fraction φ(J ). This framework became known as the μ(J )-rheology and allowed to describe the critical values μ c and φ c as particle properties that determine the jamming transition, i.e., J → 0 and φ → φ c . Note that the viscous number has initially been introduced as I v [168], which could lead to the misleading and contradicting translation of a viscous inertial number. Instead, we follow the terminology of [167] to use the symbol J for the viscous number. Even though the correlations for effective viscosities and macroscopic friction were derived from bulk averaged flow behavior of a dense suspension of neutrally buoyant particles in annular Couette-type flows, recent experimental [170] and numerical studies [21,44] have presented evidence that these frameworks also hold on the local microscale for sediment beds sheared by laminar viscous flows. In these studies, it was shown that one can infer the rheology of sheared sediment beds from depth-resolved, time-averaged profiles of total stress τ and granular pressure p s . The rheology could be inferred from these data by averaging over control volumes of only a few particle diameters thickness. This approach opens up a promising avenue to further extend the μ(J )-rheology for large-scale sediment transport modeling.
Working towards this aim will require more efforts in the future to explore the rheological behavior of sheared sediment beds for a wider range of flow conditions [171]. Tapia et al. [172] were able to show that modifying particle roughness directly impacts μ c and φ c . Aussillous et al. [173] pointed out the difficulty to measure the particle volume fraction in dense granular system, but this quantity is a key component to determine the granular stress p s in dense granular packings that enter the friction-dominated quasi-static regime [105]. Consequently, simpler empirical models that were derived from plane shearing [174] are currently employed in two-phase fluid solvers [109]. Recently, progress has been made in this regard by Vowinckel et al. [175], who studied the rheology of sediment beds sheared by pressure-driven flows by comparing the data of Aussillous et al. [173] to highly resolved data generated by means of particle-resolved simulations following the scheme of [44,49]. The presented analysis revealed good agreement of the datasets with the μ(J )-rheology in the quasi-static limit, but also highlighted the difficulty to accurately capture the rheology within the sediment transport layer.
Furthermore, it is obvious that additional effects that have been described in Sects. 3.1 and 3.2 such as cohesion, polydispersity, and non-spherical particle shapes may also have an impact on the rheology. Extending rheological correlations to incorporate these effects would be a great step forward for two-phase flow modeling. For sediment transport phenomena, it is important to recall that sediment beds sheared by a viscous fluid flow can be subdivided into two different regions [176]. At larger depths, the sediment is densely packed and the μ(J ) holds, because the dynamics of the dense suspension is governed by stresses from frictional contact. Close to the interface of the sediment bed and the clear fluid flow, i.e., the sediment transport layer, the local volume fraction decreases and the particle-particle interaction is dominated by binary collisions. Hence, Maurin et al. [177] suggested that the classic rheological frameworks may need to be extended to the more dilute regime of active sediment transport.
Another popular framework to close the system of equations that make up the two-phase flow model is the kinetic theory for granular materials. Originally, the kinetic theory was derived for heat transfer in gases, but it was shown that these principles can be extended to granular flows by defining a granular temperature as a means to quantify the kinetic energy of the sediment grains [178]. Similar to the μ(J )-rheology, the starting point was the consideration of dry granular media with low to moderate volume fractions [109]. The framework was later extended to particle transport in viscous shear flows by providing closures for the granular pressure, the collisional and frictional particle stress, the effective viscosity, the granular temperature flux, the dissipation rate, and the fluid-particle interactions [179].
Owing to the fact that the general idea of the kinetic theory framework is more suitable for low to moderate volume fractions, the original theory does not work so well for dense suspensions, where enduring contacts of the particles exist [180]. Hence, improvements of the kinetic theory are still needed to use it as a macroscopic model for sediment transport phenomena. For example, efforts have been made to extend the theory by empirically accounting for the rheology of a dense suspension [181] and to account for finite-sized particles when dealing with the stress exchange between the fluid and the particle phase. It was shown by Mathieu et al. [182] that this stress exchange term bears strong similarity to the BBO equation (12) that contains the assumptions of drag relations for isolated, spherical particles (cf. Sect. 2.1.2). These relations can be corrected to obtain consistent drag models for two-phase fluid and point particle approaches [183], but a more rigorous extension to the dense regime as it is accurately captured by the μ(J )-rheology may be needed. The limited range to which both the μ(J )-rheology and the kinetic theory are applicable, leave it up to the user to make a proper choice for the macroscopic closure depending on the physical scenario under consideration [180]. It would ultimately be desirable to extend the kinetic theory and the μ(J )-rheology in such a way that, when applied independently in a two-phase flow model, they can provide consistent results for a specific physical scenario of sediment transport.
The need for research is especially important, considering that both approaches, kinetic theory and μ(J )rheology, have already been implemented in open source CFD software packages, such as the extension of OpenFOAM [184,185] to develop a three-dimensional two-phase flow solver for sediment transport (SedFoam) [109,186]. The μ(J )-rheology has readily been used for simulations of the granular collapse in a viscous fluid [187].

Turbulence modeling for sediment transport
Another aspect that deserves further research is the issue of proper turbulence modeling. After all, the purpose of using continuum models is to investigate larger systems that can no longer resolve the smallest scales of fluid and particle motion. So far, applications involving the two-phase flow model use classical turbulence models such as the dynamic model [188] for LES [110,182] or k − and k − ω models for RANS simulations [186,187]. It has been shown, however, that the momentum balance for fluid flows over rough topography, such as aquatic bed forms, may require an extension, because the timescales of the bed morphology become much larger than the fluid timescales [189]. To this end, Nikora et al. [190][191][192] developed the double-averaging methodology (DAM) to extend the momentum balance of turbulent flows over fixed rough beds by introducing appropriate averaging operators. The averaging principal, which is similar to the Favre averaging Eq. (19) [106,108], separates the fluid from the solid phase by an indicator function. In addition, the DAM rigorously decomposes a given quantity into a value averaged in space and time and a temporal and a spatial fluctuation. This decomposition gives rise to additional stress terms that originate from the flow over the heterogeneous terrain of the sediment bed, i.e., bed forms.
The framework was later enhanced for turbulent flows over mobile granular sediment beds [193]. Subsequently, the importance of the different components entering the extended momentum balance was illustrated by means of particle-resolved DNS of mobile sediment transported in a turbulent open-channel flow [54,55]. The DAM has recently been extended from the momentum balance to the second-order velocity moments [194], which extends classical RANS modeling from the kinetic energy budget to the turbulent double-mean, form-induced, and spatially averaged energy budgets. Using the same data set of Vowinckel et al. [54,55], an analysis of all relevant terms entering the extended kinetic energy budget has been presented by Papadopoulos et al. [195]. This analysis showed that, indeed, the extended energy budget needs to be considered to close the budget of the entire system. However, the DAM inherits the closure problem from the classical RANS equations, which introduces additional stresses arising from the averaging of the nonlinear Navier-Stokes equations. To eventually apply the DAM for two-phase flow models as a new type of double-averaged Navier-Stokes (DANS) equations, closures will be needed to obtain unique solutions for the extended system of equations.

Machine learning
Machine learning (ML) has been used since the 1950s and was originally developed for gaming and pattern recognition [196]. Generally speaking, this technique is meant to translate a given input vector into a desired output without explicitly knowing the connection between the two. The algorithm for that translation, hence, remains a "black box" where the output is an approximate solution that is not necessarily based on exact physical laws. There are different strategies depending on whether there exists a labeled data set of expected results. While unsupervised learning tries to make sense of unlabeled data by extracting features and/or patterns, an algorithm that compares its results to a labeled data set is known as supervised learning. A key feature of supervised ML is that the algorithm can learn from training data, by constantly improving its accuracy as the Fig. 11 Sketch of a deep neural network to map a given input I to a desired output O with corresponding weights w, biases b, and activation functions a, where the superscript indicates the index of the hidden layers and the subscripts indicates the node in the present and preceding layer, respectively. The total number of degrees of freedom is determined by the number of hidden layers i and the number of nodes of a given layer, here j, k, and l (modified after Ling et al. [203]) computed output is compared to the expected results of the labeled data set. The observed deviation can then be translated into a correction algorithm that optimizes the free parameters of the black box model to minimize the error. The accuracy of supervised ML algorithms should, therefore, increase with the amount of training data fed into the optimization procedure.
One popular class of black box models used for supervised learning is a network of interconnected layers of so-called artificial neurons [197], a technique that became known as Artificial Neural Networks (ANN). An example of an ANN is given in Fig. 11, where the part that is considered as a "black box" is marked in the light yellow color. For ANNs, the input, denoted as I, is passed as a weighted sum through a network of interconnected neurons that all contain weights w i l, j , biases b i l , and nonlinear activation functions a i l . In case the number of hidden layers is larger than one, i.e., i > 1, the network is called a deep neural network (DNN). The weights and biases all represent degrees of freedom that can be optimized using iterative procedures to obtain the desired output. After a successful training phase, the optimized algorithm needs to be tested against a validation data set that has not yet been used for training. This step is needed to assess the predictive capacity of the ML algorithm.
Since the early 1990s, ANNs have been applied for modeling techniques in the field of hydrology and environmental engineering including runoff and flood forecast from a watershed, water quality, wastewater treatment, ground water flows and, in particular, the discharge of suspended sediment [198][199][200][201][202]. This strategy was especially useful at the scale of entire watersheds, because the underlying physical principles of the governing processes could neither be accessed by measurements nor resolved by the modeling techniques at the time and it will be unlikely to access such scales with the numerical techniques outlined in Sect. 2 in the foreseeable future.
Nevertheless, exploiting ML techniques for fluid dynamics applications has become increasingly popular because of its potential to connect the different scales typically encountered in environmental fluid dynamics. For example, turbidity currents can range from particle-particle interactions on the scale of micrometer and milliseconds to their total runout time and length at the order of several days and O(1000km), respectively [2]. Bridging this disparity in scales will inevitably require heuristic rules/closures [204], so that the unknown physical connection between the scales can in parts be replaced with modern ML techniques. Owing to the increasing computing capacity, seminal work has been performed in recent years to advance subgrid-scale turbulence modeling from highly resolved DNS to larger-scale RANS simulations of turbulent flows [203,205,206]. In the context of multiphase flows, further efforts have been made to derive subgrid-scale models for two-phase flow models of bubbly flows [207] and the percolation of fine-grained sediments into a porous sediment bed [208].
Even though the connection between the scales may not always be explicitly known, it is important to highlight that fluid dynamics principles of the multiscale phenomena such as conservation of mass and momentum cannot be ignored [209]. Otherwise, the results provided by ML algorithms may not hold for validation, and a physical interpretation is no longer possible [210]. This is especially true for highly resolved simulations of sediment transport phenomena such as particle-resolved DNS. Consequently, these constraints  [212]). The flow is driven by a constant pressure gradient through an array of fixed spheres. Iso-surfaces represent a constant value of streamwise velocity may render DNN approaches in parts inadequate due to the large number of free parameters that bear the risk of overfitting results. That is, the accuracy of the validation phase decreases with an increase of free parameters [211]. On a similar note, it is important to realize that, strictly speaking, ML algorithms provide efficient interpolation for the range defined by the training data. Hence, extrapolating beyond the bounds set by these data may yield inaccurate results that can no longer be trusted or interpreted in a physical sense [204].
Consequently, the most important task for the use of ML in the context of sediment transport applications will be to properly reduce the vast amount of data from highly resolved numerical simulations to a meaningful set of input parameters that can be handled by the learning algorithm [213]. This will be a key for its efficient use, because it is known that the complexity of the problem increases exponentially with the number of input parameters. This dilemma is referred to as the curse of dimensionality [214]. Two very promising frameworks that derive hydrodynamic forces acting on individual particles in a dense suspension have recently been proposed as the pairwise interaction extended point-particle (PIEP) [212] and the microstructure-informed probability-driven point-particle (MPP) model [18], respectively. Both models have in common that data from particle-resolved DNS of a flow through an array of fixed spheres is used to derive the effect of particle hiding and shading on the resulting fluid drag acting on individual particles in that same configuration (Fig. 12). To this end, the dimensionality is reduced by superposing the effect of pairwise particle interaction, so that the overall drag correction of a single particle becomes the sum of the effects experienced from the neighboring particles [69]. The required input of the ML algorithm then becomes the Reynolds number based on the mean flow, the particle diameter, and the local value of the volume fraction as well as its position and its nearest neighbors to obtain a given drag force. The algorithm is optimized by using the drag force computed by the particle-resolved DNS as training data to quantify the pairwise interaction. This interaction can be mapped out in terms of force maps [212] or probabilities of finding a neighbor in a certain distance for a given recorded drag force [18]. Nonlinear least square methods, such as the Levenberg-Marquardt algorithm, are then used to optimize a set of open parameters that enter, e.g., a spherical harmonic expansion, to recover the hydrodynamic forces obtained from particle-resolved DNS with satisfactory agreement. Subsequently, drag forces in point-particle EL simulations can be computed for a given particle arrangement with high fidelity.
The huge potential of this approach can be further enhanced towards more physical realism regarding the situation of sediment transport in naturally occurring flows. For example, additional campaigns of these types of investigations could involve ensembles of polydisperse particles, allow for particle rotation, include effects of fluid shear, particle acceleration, etc. More broadly speaking, these tools can also be used to incorporate the effects of cohesive and non-spherical particles described in Sects. 3.1 and 3.2 into continuum-type frameworks discussed in Sect. 3.3. Nevertheless, it is also important to keep in mind the risks of the excessive use of ML techniques. Those are (i) the curse of dimensionality, (ii) overfitting the ML algorithm, (iii) extrapolating the data beyond the bounds of the training data, and (iv) the loss of physical reasoning for a black-box model that has little or no input from physics-informed processes. Hence, while ML has the potential to serve as a tool that helps to bridge the gap in between scales of different simulation approaches, care should be taken as to how robust and reliable those newly developed tools are.

Concluding remarks
Much progress has been made in sediment transport modeling over the last decades. The recent advancements in computational fluid dynamics coupled with suitable discrete element methods (CFD-DEM) already allow to investigate a host of sediment transport phenomena on different scales. These achievements are strongly coupled to the increase of computational power. It is now possible to perform particle resolved simulations of turbulent open-channel flows that resolve the fluid motion on the grain scale and thereby reproduce bed forms such as ripples and dunes with high fidelity. These highly resolved simulations can readily be utilized to understand the governing processes of bed form evolution. Nevertheless, it will be important to further enhance the CFD-DEM technique towards a better representation of natural environments. Two key features that address this issue have been discussed in this review. These are the impact of inter-particle cohesion and particle shape on the sediment mobility.
It is equally important, however, to translate these processes into constitutive equations that feed into larger-scale models that are computationally more efficient. Whenever possible, it is desirable to derive these constitutive equations by means of physically based scaling laws. A prime example discussed in this review is the rheological behavior of dense particle systems sheared by viscous flows. Promising rheological scaling laws have already been proposed and successfully implemented into two-phase fluid models that treat the sediment phase as a continuum. To date, different approaches exist for the rheology of sediment beds that have been derived for either more dense or rather dilute particle systems, respectively, but have limitations regarding certain extreme flow regimes. Hence, more work is needed to arrive at a consistent set of tools to increase the predictive capacity of two-phase fluid models. In addition, more work will be needed to incorporate the effects of the long-term morphological evolution that may exceed the timescales considered in classical two phase fluid models.
Owing to the complex interactions in a fluid-particle mixture, a physically based derivation of the desired scaling laws may not always be possible in a straightforward manner. Physics-informed machine learning techniques may, therefore, become a more important alternative to derive these closures from highly resolved data, but care should be taken regarding its limitations to retain the physical reasoning of data-driven machine learning algorithms. Nevertheless, combining these three efforts from particle-resolved simulations, enhanced upscaling strategies, and comprehensive two-phase fluid models has the potential to not only advance our understanding of sediment transport from a scientific point of view, but also become an established tool for engineering-type practices in the future.