Interactive comment on “ SedFoam-2 . 0 : a 3 D two-phase flow numerical model for sediment transport ” by Julien Chauchat et al

This paper presents a 3D two-phase flow numerical model for sediment transport (SedFoam-2.0) in detail, including the mathematical formulation and the numerical implementations. The authors newly include the mixing length turbulence model, the k − ω model, and dense granular flow rheology into SedFoam. The main purpose is to provide a comprehensive numerical framework that solves the two-phase flow equations in three dimensions with the capability to select different combinations of turbulent model and granular stress model for sediment transport. This paper is well written and pleasant to read. The reviewer suggests acceptance after minor revisions.


Introduction
Sediment transport is the main process that drives the morphological evolution of fluvial and coastal environments.Consequently, the ability to predict sediment transport is a major societal issue for the management of natural sys-tems in order to limit and prevent the impacts related to extreme events exacerbated by climate change and human activities such as the construction of hard structures (dams, harbors, dikes, etc.), land reclamation, and dredging.Addressing these issues requires the development of comprehensive models that account for the variety of complex hydrosedimentary processes such as particle interactions with hydrodynamic and flow turbulence or particle-particle interactions due to collisions or frictions.However, these complex phenomena are only poorly understood at present and they are incompletely integrated into engineering tools to predict the coastal and river morphodynamics.As a result, our prediction performance is limited.Improving these models is urgently needed for land settlement decision-makers for the management of water resources and environmental issues.
The processes at work in sediment transport are numerous.In classical sediment transport definitions, particles can be transported as suspended load, i.e., without contact with the sediment bed, or as bed load, i.e., with permanent or intermittent contact with the sediment bed by rolling, sliding, or saltation (Fredsoe and Deigaard, 1992).In the suspended load, sediment concentration is quite low and the sediment suspension is driven by their interactions with turbulent eddies.In the near-bed region, the sediment concentration increases drastically and can reach values as high as ∼ 60 % in volume fraction.As a result, the particle-particle interactions such as collisions or frictional contacts become dominant.When the granular particles are agitated by a strong shear stress, a thick layer of particles are mobilized and transported Published by Copernicus Publications on behalf of the European Geosciences Union.
above the static bed.This regime, the so-called sheet-flow regime, can significantly contribute to the sediment transport and the morphological evolution of rivers (for example, under extreme flood conditions, Hanes and Bowen, 1985) and in the coastal zone (especially in the surf and swash zones, e.g., Lanckriet et al., 2014;Aagaard et al., 2002).
From the modeling perspective, the classical modeling approach consists of dividing the physical domain into two sublayers.The upper layer corresponds to the water column in which depth-integrated or depth-resolving Reynoldsaveraged Navier-Stokes equations are solved and the sediment concentration is assumed to be dilute, in which the sediment particles are treated as passive scalar with a settling velocity difference with the fluid phase.The lower nearbottom bed-load layer is solved by using a bed shear-stressbased empirical formula for the sediment flux (e.g., Meyer-Peter and Muller, 1948) coupled with the sediment mass conservation equation (Exner equation, Exner, 1925).The latter allows us to predict the bed morphological evolution.These two layers are dynamically coupled using empirical sediment vertical fluxes: deposition and pick-up fluxes (Van Rijn, 1984), which are also parameterized based on the bed shear stress.The single-phase model has been widely used because of its simplicity and computational efficiency, and it has been integrated into meso-to large-scale models such as the Delft3D (Lesser et al., 2004;Hu et al., 2009), XBeach (Roelvink et al., 2009), and ROMS (Warner et al., 2008).Despite their ability to predict long-term regional morphology in littoral zones, there are several challenges to model sediment transport using such a single-phase methodology.Firstly, major transport in sheet flow occurs within about 20-50 grain diameters above the bed.Hence, the resolved suspended-load layer only accounts for a minor portion of transport while the majority of the predicted transport relies on the empirical bed-load parameterization.Secondly, most of the existing formula for bed-load transport and suspension flux are developed for steady flow and hence their applicability to the highly unsteady wave-induced transport is questionable (Yu et al., 2012).Thirdly, the particle-particle collisions and frictions in the high concentration regions are assumed to be negligible, and entrained sediments are assumed to be suspended immediately by the flow turbulence.Fourthly, the bed erosion/deposition and its impact on the sediment concentration distribution and the carrier flow fields are largely neglected.In addition, simple parameterization solely on the bed shear stress may be insufficient; for example, the role of pressure gradient in sediment transport has been identified for extreme events such as storms (Foster et al., 2006;Cheng et al., 2017a).Addressing these issues requires the development of comprehensive models that account for the variety of complex hydrodynamics and sediment transport processes on a regional-scale setting (e.g., Lesser et al., 2004;Roelvink et al., 2009).Because sediment transport occurs very close to the bed, effective parameterizations of sediment transport are needed.In the past decade, significant progress has been demonstrated to utilize detailed numerical models to understand sheet-flow processes, and effective parameterizations for coastal sediment transport have been developed (van der A et al., 2013).To further tackle more complex sediment transport problems, such as bedform evolution, scour, bank erosion, and dune erosion, further expansion of these models through a community effort is urgently needed.The development of more comprehensive sediment transport models integrating the complexity of the underlying physical coupling mechanisms is the main goal of the open-source community model presented herein.
During the past two decades, an increasing amount of research efforts have been devoted to develop two-phase flow models for sediment transport (see a brief summary in Table 1).In this two-phase flow approach, dynamical equations are solved for both the fluid phase (water) and the particle phase (sediment), with the latter being seen as a continuous phase dispersed in the fluid.The two-phase flow approach gives a general modeling framework that potentially allows us to take into account almost all the physical processes involved in sediment transport, such as fluidparticle interactions, turbulence modulation, and particleparticle interactions.From a technical point of view, the difficulty is in solving two non-linearly coupled Navier-Stokes type of partial differential equations.From a theoretical perspective, it is difficult to incorporate various physical processes that take place on smaller scales than the averaging scale, which is used to derive the two-phase equations.With closures for particle-particle interactions, flow turbulence, and turbulence-sediment interactions, the Eulerian two-phase model excels the single-phase sediment transport models in several aspects, and provides us with new insights into sediment transport mechanisms.
The key closures in the two-phase flow sediment transport models are flow turbulence and granular stress closures.In terms of the turbulence closure, the first one that has been tested was a mixing length model by Jenkins and Hanes (1998) followed by others (e.g., Dong and Zhang, 1999;Revil-Baudard and Chauchat, 2013).The other turbulence model that has been used is the k − ε model (e.g., Hsu et al., 2003;Longo, 2005;Bakhtyar et al., 2009).The k − ω model has been tested by Jha and Bombardelli (2009) and Amoudry (2014), and the Reynolds stress model has been tested by Jha and Bombardelli (2010).Concerning the granular stress models, the first model that has been tested by Hanes and Bowen (1985) is the empirical rheology of Bagnold (1954).This rheology has then been used by Dong and Zhang (1999) and Bakhtyar et al. (2009) for oscillatory sheet-flow applications.Jenkins and Hanes (1998) were the first to apply the kinetic theory of dense granular flows in a two-phase flow model for sediment transport.The kinetic theory has been used quite extensively by different authors to study sediment transport (e.g., Hsu and Liu, 2004;Berzi, 2011;Berzi and Fraccarollo, 2013;Cheng et al., 2017a).In this approach, the particle stress associated with particle-particle collisions

Authors
Turbulence model Particle stress Asano (1990), Li and Sawamoto (1995), and Dong and Zhang (2002) mixing length Bagnold Jenkins and Hanes (1998) mixing length kinetic theory Revil-Baudard and Chauchat (2013) mixing length granular rheology Li et al. (2008) k − L Bagnold Bakhtyar et al. (2009) k − ε Bagnold Hsu and Liu (2004), Chauchat and Guillou (2008), and Amoudry et al. (2008) k − ε kinetic theory Yu et al. (2010) and Cheng et al. (2017a) Amoudry (2014) and Jha andBombardelli (2009, 2010) k − ω kinetic theory Lee et al. (2016) k − ε granular rheology is modeled by the fluctuation energy of the particle phase (or granular temperature).Various models were developed to model the granular temperature.Jenkins and Hanes (1998) first applied kinetic theory for dry granular flows to sheet flow, and the granular temperature transport equation was later extended to consider the fluid-sediment turbulence interactions (e.g., Hsu and Hanes, 2004;Chauchat and Guillou, 2008).The aforementioned kinetic theory considered dense collisions of particles, while the streaming effect of particle random motions was missing in the dilute concentration regime.A further extension to include the streaming effects was developed by Lun and Savage (1987) and Ding and Gidaspow (1990), and it was implemented into a more complete two-phase model by Cheng et al. (2017b).In contrast to solving the transport of granular temperature, Jha and Bombardelli (2010) used a mixing length concept for the particle phase, and a simpler algebraic model for the granular temperature was used with success.More recently, the dense granular flow rheology initially proposed by GDRmidi (2004) for dry granular flows has been used for sediment transport applications in the laminar flow regime by Ouriemi et al. (2009) and later to turbulent flow conditions by Revil-Baudard and Chauchat (2013), Chiodi et al. (2014), and Lee et al. (2016).
Due to the complexity of the model formulation, most of the existing two-phase models are based on the Reynoldsaveraged approach and simplified into one-dimensional form (e.g., Hanes and Bowen, 1985;Jenkins and Hanes, 1998;Hsu et al., 2003;Revil-Baudard and Chauchat, 2013), with a few exceptions (two-dimensional models; Chauchat and Guillou, 2008;Bakhtyar et al., 2009;Amoudry and Liu, 2009).Only very recently, a three-dimensional (3-D) large-eddy simulation two-phase flow model has been applied for sheet-flow sediment transport (Cheng et al., 2017b).This 3-D numerical model is based on the open-source CFD toolbox Open-FOAM.Cheng et al. (2017a) have applied SedFoam using the kinetic theory of granular flows and the k − turbulence model to reproduce oscillatory sheet flows of fine, medium, and coarse sand (O'Donoghue and Wright, 2004).The purpose of the present contribution is to follow up on the work of Cheng et al. (2017a) by adding new capabilities to the open-source model SedFoam.In particular, the mixing length turbulence model and dense granular flow rheology used by Revil-Baudard and Chauchat (2013) and Chauchat (2017) for sheet flows have been implemented.In addition, we implemented and tested the k − ω turbulence model for twophase flow sediment transport modeling purposes.Our final goal is to provide a comprehensive numerical framework that solves the two-phase flow equations in three dimensions with the capability to select different combinations of turbulence and granular stress models for sediment transport applications.By disseminating the numerical model in the open-source framework, in the long run, we expect new capabilities will be added to the model by the scientific community.We strongly believe that developing such an opensource community model is the only effective way to make significant progress.
The focus of the present work is to demonstrate the capabilities included in SedFoam-2.0 to model sediment transport.In particular, the comparison of different combinations of granular stress and turbulence models in the same numerical framework is presented here for the first time.
The paper is organized as follows, in Sect. 2 the mathematical formulation of the model is presented together with the closures for drag, turbulence model, and granular stress models.In Sect.3, the semi-discretized form of the equations and the velocity-pressure coupling algorithm are presented in detail.In Sect.4, four test cases are presented: three of them are one-dimensional vertical problems namely sedimentation, laminar bed load, and sheet flow for which experimental data or analytical solutions exist.The last test case on the scour at an apron is used to illustrate the multidimensional capabilities of SedFoam-2.0.Finally, summaries and conclusions are drawn in Sect. 5.

Mathematical formulation
The mathematical formulation of the Eulerian two-phase flow model is obtained by averaging local and instantaneous mass and momentum conservation equations over fluid and dispersed particles.Different averaging operators can be www.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017 used, ensemble averaging (Drew, 1983) or spatial averaging (Jackson, 2000), and provided that the mathematical derivation is done properly, the different approaches should lead to the same conservation equations (Zhang and Prosperetti, 1997;Jackson, 1997).The resulting governing equations can be considered as the counterpart of the clear fluid Navier-Stokes equations for single-phase flow.In order to apply these equations to turbulent flow, in which turbulent motions are generated by flow shear much larger than the grain scale, additional turbulence averaging or filtering is required.In the present model, turbulence-averaged Eulerian two-phase flow equations are derived by following a similar procedure presented in Hsu et al. (2003) and Hsu and Liu (2004).

Turbulence-averaged two-phase flow governing equations
The mass conservation equations for the particle phase and fluid phase are written as the following: where α and β (β = 1 − α) are the particle and fluid volume concentrations, u a i , u b i are the sediment and fluid-phase velocities, and i = 1, 2, 3 represents streamwise, spanwise, and vertical components, respectively.The momentum equations for fluid and particle phases can be written as the following: where ρ a , ρ b are particle and fluid density, respectively, g i is the gravitational acceleration, f i is an external volume force, and p is the fluid pressure.f i is the external force that drives the flow.The fluid stress τ b ij includes fluid grain-scale (viscous) stress and fluid Reynolds stresses (see Sect. 2.2), and p a , τ a ij are particle normal stress and shear stress (see Sect. 2.3).The last two terms on the right-hand side (RHS) of Eqs.(3) and (4) are momentum coupling between the fluid phase and particle phase through drag force, where K is the drag parameter.In particular, the second to last term represents averaged drag force due to mean relative velocity between fluid and particle phases, while the last term represents the fluid turbulent suspension term, also called drift velocity by Simonin (1991).This term is due to the correlation of sediment concentration and fluid velocity fluctuations, and the gradient transport assumption is adopted here for its closure.Hence, ν b t is the turbulent viscosity to be calculated using a turbulence closure, and S US = 1/σ c is inverse of the Schmidt number.This term is equivalent to the turbulent suspension flux of the Rouse profile in the two-phase flow formalism (see Chauchat, 2017, Appendix 1 for a detailed demonstration).Other forces such as the lift force or the added mass force could play a role in sediment transport: according to Jha and Bombardelli (2010), the lift force in dilute suspended sediment transport only represents 4 % of the drag force and the added mass force can be on the order of 10 % in the near-bed region.The influence of the added mass force would require further investigation that is beyond the scope of the present paper.
The drag parameter K is modeled following Schiller and Naumann (1933): where d eff = ψd is the effective sediment diameter, in which ψ is the shape factor and d is the particle diameter.The hindrance function β −h Exp represents the drag increase when the particle volume concentration increases.h Exp is the hindrance exponent that depends on the particulate Reynolds number (Di Felice, 1994).For simplicity, the value of h Exp is assumed to be constant (default value is 2.65), and its value can be specified from the constant/transportProperties file in SedFoam-2.0.This hypothesis is valid for particulate Reynolds numbers lower than unity or larger than 300; within this range the exponent value decreases down to h Exp ≈ 2. The drag coefficient C d is calculated as , where ν b stands for the fluid kinematic viscosity.This drag model can be chosen by the keyword "GidaspowSchillerNaumann" in the file constant/interfacialProperties and it is especially well adapted for dealing with suspended particles.For situations in which the fluid flow inside the porous sediment bed has to be accurately solved, other drag models are available in SedFoam-2.0but as they are not relevant to the test cases investigated in this paper, their description is omitted.The major issue in developing a Eulerian two-phase flow model is to provide closure laws for turbulence closures and granular stress models.This will be extensively discussed in the following subsections, where particularly different modeling options available in SedFoam-2.0are presented.

Fluid-phase shear stress
Because the present model equations are obtained by averaging over turbulence, the fluid stresses consist of a large-scale component R bt ij (i.e., Reynolds stress) and a grain-scale stress r b ij , which includes the viscous stress and an additional effect due to fluid-particle interaction on the grain scale.The total fluid stress is written as where ν b Eff = ν b t + ν mix is the fluid-phase effective viscosity with ν b t being the eddy viscosity, and ν mix is the mixture viscosity.S b ij is the deviatoric part of the fluid-phase strain rate tensor defined as The Reynolds stress tensor R bt ij is modeled as and the viscous stress tensor is modeled as In SedFoam-2.0,several different viscosity or turbulence closures are implemented, and these models can be selected according to specific flow conditions ranging from laminar to turbulent flows, and in particular, the mixture viscosity can be selected in combination with the granular rheology model for the granular stresses (see Sect. 2.3.2).For the turbulent eddy viscosity, modified turbulence closures are implemented for sediment transport applications, and they are presented in Sect.2.2.2.

Mixture viscosity
The mixture viscosity model mostly depends on the particlephase volume concentration.Four different models are available in SedFoam-2.0.In the pure fluid model, the mixture viscosity is equal to the fluid one: ν mix = ν b .This model is selected by default when using the kinetic theory of granular flows and can be selected by setting none for the keyword FluidViscosityModel in the file constant/granularRheologyProperties.
The Einstein model (Einstein, 1906) is valid for very dilute situations (α < 0.01): The phenomenological model proposed by Krieger and Dougherty (1959) is valid for very dense situations: where α max is the maximum volume concentration and n is an empirical exponent usually taken as n = 2.5α max for consistency with Einstein's model at low volume concentration.
The model proposed by Boyer et al. (2011) is also empirical but has been obtained based on detailed rheological experiment.It is consistent with Einstein's model for low volume concentration and Krieger-Dougherty's model at high volume concentrations: The choice of mixture viscosity model is made in the file constant/granularRheologyProperties through the keyword FluidViscosityModel and is only available when the granular rheology is chosen.

Turbulence modeling
As discussed above, the turbulence-averaged formulation requires a closure for the eddy viscosity.Three turbulence models are available in SedFoam-2.0:a mixing length model (only valid for 1-D configuration), the k − ε model from Cheng and Hsu (2014) and Cheng et al. (2017a), and a k − ω turbulence model introduced in the present contribution.

Laminar
For laminar-flow applications, the turbulence model is turned off by setting ν b t = 0 and k = 0; however, the mixture viscosity model can be selected to account for the sediment effect on the mixture viscosity; thus the effective fluid viscosity is calculated as ν b Eff = ν mix .

Mixing length (1-D only)
In the mixing length approach, the eddy viscosity is modeled using a simple algebraic equation: where κ is the von Karman constant and the exponent 1.66 has been proposed by Chauchat (2017) based on matching experimental data from Revil-Baudard et al. (2015).This model is only working in the 1-D configuration for which the direction of gravity is y.This turbulence model has been implemented mostly for compatibility with earlier works.

Keyword name Keyword value
FluidViscosityModel none Einstein KriegerDougherty BoyerEtAl Eq. ( 11) Eq. ( 12) Eq. ( 13) the turbulent eddy viscosity ν b t is calculated by where C µ is an empirical coefficient (see Table 3).The turbulent kinetic energy (TKE) k is computed from the solution of Eq. ( 17), appropriate for sand particles in water: The above k equation is similar to the clear fluid k−ε closure, except that the last two terms on the RHS in Eq. ( 17) take account of the sediment damping effect on the carrier-flow turbulence through drag (the fourth term) and density stratification (the last term).In the drag-induced damping term, the parameter t mf is introduced to characterize the degree of correlation between particles and fluid velocity fluctuations and it can be quantified by the Stokes number St (Benavides and van Wachem, 2008): where t p = ρ a /(β K) is the particle response time and t l = k/(6ε) is the characteristic timescale of energetic eddies.Danon et al. (1977) and Chen and Wood (1985) proposed an exponential function for t mf , which is also used in Kranenburg et al. ( 2014) and Cheng et al. (2017a): where B is an empirical coefficient.The last term in the TKE Eq. ( 17) represents the buoyancy term.For typical sediment concentration with an upward decaying profile, this term represents the well-known sediment-induced stable density stratification that provides another source of turbulence attenuation.
Finally, the balance equation for the rate of turbulent kinetic energy dissipation ε is written as As discussed in Hsu et al. (2004), due to a lack of comprehensive experimental data, the coefficients associated with the present two-equation closure are adopted from their clear fluid counterparts.The coefficient C 3ε in the ε Eq. ( 20) is chosen to be 1.2.For the coefficient associated with the buoyancy term, C 4ε = 0 is used for the stably stratified condition, while it is set to 1 for the unstably stratified condition.Table 3 summarizes the model coefficients.These coefficients were shown to work well for typical medium to coarse sand transport (Hsu et al., 2004;Yu et al., 2010;Cheng et al., 2017a).Furthermore, it was found that the coefficient B (see Eq. 19) is sensitive to the model result and thus the coefficient B is chosen to be a free parameter to be calibrated with measured data.

k − ω model
It is well-known that the original k − model has been derived for high Reynolds number flows and is not very accurate to describe transitional flows such as the situation of the flow reversal in a wave boundary layer (Guizien et al., 2003).For this situation and for near-wall treatment, the k−ω model is more suitable and more stable than the k − model (Guizien et al., 2003).Another physical situation in which a k − ω model works better than a k − model is in the presence of an adverse pressure gradient such as the downward facing step or at the upstream side of an obstacle (Menter, 1994;Wilcox, 2008).In order to test the influence of the turbulence model, a two-phase k − ω model is introduced in the present contribution, which is very similar to those of Jha and Bombardelli (2009) and Amoudry (2014).
The turbulent eddy viscosity ν b t is calculated as Following the same method as for the two-phase k − ε turbulence model, the fluid TKE and the specific turbulent energy and the equation for ω reads as The different coefficient values can be found in Table 4. Similar to the two-phase k − ε model, the coefficients are equal to the clear fluid model.Only C 3ω and C 4ω need to be determined.According to the numerical experiments described in Sect.4, the coefficient C 3ω is chosen to be equal to 0.35 and C 4ω = 0 is used in stably stratified situations, while it is set to 1 for unstably stratified situations.Due to different model configurations in k − ε and k − ω models, the coefficient B may be different.However, a similar sensitivity in the B coefficient (see Eq. 19) to the model results is observed, and it is left as the only free model calibration parameter.
The turbulence model can be selected using the RASModel keywords in the file constant/RASProperties and specific parameters of the two-phase turbulence models are set in the file constant/twophaseRASProperties (see Table 5).
Note that Tables 3 and 4 present the default coefficient values used in our implementation, which has been validated extensively with measured sediment transport data.However, they can be changed by setting the keyword and values in the files constant/RASProperties and constant/twophaseRASProperties (see notations list in Appendix A).

Particle-phase stress
In sediment transport applications, the particle stresses are important mechanisms to support a particle's immersed weight in concentrated regions of sediment transport (Hsu et al., 2004;Cheng et al., 2017a).In these regions, the momentum exchanges due to particle collisions and/or enduring contacts exert dispersive stresses on a collection of particles.The particle-phase stress tensor can be split into the normal and off-diagonal components that correspond to the particle pressure p a and the particle shear stress τ a ij .Although details vary with the model selection, the particle normal stresses (or pressure) can be generally classified into two contributions: a shear-induced or collisional component (super-script "a") and a permanent contact component (super-script "ff ") (Johnson and Jackson, 1987): The first term of the particle pressure is due to enduring contact in the highly concentrated region, where sediment bed is quasi-static or immobile.This normal pressure increases rapidly when the sediment concentration is close to its maximum packing limit and prevents unphysical sediment concentration in the sediment bed.Thus, this element is important to model a full transport profile including the quasi-static sediment bed.The permanent contact component p ff is calculated as where α Fric min = 0.57, α max = 0.635 for spheres and Fr, η 0 , and η 1 are empirical coefficients.Following Cheng et al. (2017a), the values are set to Fr = 0.05, η 0 = 3, and η 1 = 5.Again, these coefficients can be reset in file constant/kineticTheoryProperties or constant/granularRheologyProperties by specifying the keywords and values (see notations list in Appendix A).
In the modern sediment transport modeling framework, two major threads of modeling approach for shearinduced/collisional particle normal stress and shear stress are kinetic theory of granular flows and dense granular flow rheology.They are implemented in this version of SedFoam-2.0(see Sect. 2.3.1 and 2.3.2).As a result of different closures for particle normal stresses p a , the closure for the total particle shear stress τ a ij also varies, and they are described in the next two sections.

Kinetic theory of granular flows
In the kinetic theory model, intergranular interactions are assumed to be dominated by binary collisions for low to moderate sediment concentration, and the collisional shear stresses www.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017

Keyword name Keyword value
RASModel laminar twophaseMixingLength twophasekEpsilon twophasekOmega (see Table 3) (see Table 4) are quantified by particle velocity fluctuations represented by the granular temperature .The model is originally developed for dry granular flow and consists of smooth, slightly inelastic, spherical particles (Jenkins and Savage, 1983;Lun and Savage, 1987;Lun, 1991).Here, we adopt the model suggested by Ding and Gidaspow (1990), which takes into account the fluid phase.The balance equation for granular temperature is written as where the first term on the RHS is the production of granular temperature, q j is the flux of granular temperature, γ is the energy dissipation rate due to inelastic collision, and J int is the production (or dissipation) due to the interaction with the carrier fluid phase.
In the 1980s, dense-phase kinetic theory of gases (Chapman and Cowling, 1970) was applied to granular flow by many researchers (Chepurniy, 1984;Jenkins and Savage, 1983;Savage, 1988).Here, we adopt the closure of particle pressure proposed by Ding and Gidaspow (1990): where e is the coefficient of restitution during the collision.
With the binary collision assumption adopted in the kinetic theory of granular flow, the radial distribution function g s0 is introduced to describe the crowdedness of a particle.In this study, we use the radial distribution function for dense rigid spherical particle gases of Carnahan and Starling (1969): Following Gidaspow (1994), the particle collisional stress is calculated as where S a ij is the deviatoric part of sediment-phase strain rate tensor: Through the kinetic theory, the particle shear viscosity is calculated as a function of granular temperature and radial distribution function: Similarly, the bulk viscosity is calculated as The closure of granular temperature flux is assumed to be analogous to Fourier's law of conduction: where D is the conductivity of granular temperature, calculated as The dissipation rate due to inelastic collisions is calculated based on that proposed by Ding and Gidaspow (1990): Due to the presence of the carrier fluid phase, carrier-flow turbulence can also induce particle fluctuations.Following Hsu et al. (2004), the fluid-particle interaction term can be expressed as In order to extend the model capability to resolve the quasi-static or immobile sediment bed, the shear stress due to frictional contact is modeled as where ν a Fr is the frictional viscosity.By following Srivastava and Sundaresan (2003), which combined the frictional normal stress from Johnson and Jackson's model (Eq.25) and the frictional viscosity from the Schaeffer (1987) model, the friction viscosity is calculated by where a constant friction angle θ f is used.This frictional shear viscosity model has the capability to capture the transition from solid-like behavior to fluid-like behavior of the sediment bed.A small number of D small = 10 −10 s −1 is added in the denominator to ensure numerical stability when shear rate in the sediment bed becomes zero.In sediment transport, the frictional component of particle pressure and particle shear stress play a definite role to ensure the existence of an immobile sediment bed and a low mobility layer of enduring contact can be modeled (Hsu et al., 2004).The total shear stress τ a ij can be calculated as a sum of the collisional-kinetic component (τ a ij ) and a frictional component (τ ff ij ): All the simulations presented in this paper have been obtained using the closures summarized in Table 6.Other models are available in OpenFOAM and it would be very easy to implement new ones.

Dense granular-flow rheology
The other alternative for modeling the particle-phase stress proposed in SedFoam-2.0consists of the dense granular-flow rheology or the so-called µ(I ) rheology (GDRmidi, 2004;Forterre and Pouliquen, 2008).Such an approach has been used with some success by Revil-Baudard and Chauchat (2013) and Chauchat (2017) to model turbulent sheet flows.Contrary to the kinetic theory of granular flows, the dense granular-flow rheology is phenomenological, and it is based on dimensional analysis.Instead of separating the collisional shear stress and frictional shear stress, the total particle-phase shear stress is related to total particle pressure p a by a dynamic friction coefficient µ (Jop et al., 2006): The dynamic friction coefficient µ depends on the dimensionless controlling number "I ".Depending on the local Stokes and particulate Reynolds numbers, the regime of the granular-flow rheology can change from free fall or grain inertia regime to viscous and turbulent regimes (Andreotti et al., 2013), and the definition of the controlling parameter "I " is different accordingly.In SedFoam-2.0, the grain inertia and viscous regimes have been implemented and can be selected using the keywords "MuI" and "MuIv", respectively.To be consistent with the definitions used in the kinetic theory, we still introduce the particle shear viscosity µ a and frictional shear viscosity ν a Fr .However, we simply set µ a = 0, and then define the frictional shear viscosity alone as (Chauchat and Médale, 2014) where S a is the norm of the shear rate tensor (Eq.30).D small is the regularization parameter, which is introduced to avoid singularity.It is set so that D small = 10 −6 s −1 for all the simulations of granular rheology presented herein except stated otherwise.

Viscous regime
In the viscous regime, the friction coefficient µ depends on the viscous number I v = ∇u a ν b /(ρ b p a ) and is calculated as For neutrally buoyant beads, the typical values are as follows: µ s = 0.32, µ 2 = 0.7, and I 0 = 0.005 (Boyer et al., 2011).The viscous regime occurs when the Stokes number, St = d √ ρ a p a /(ρ b ν b ) as defined by Cassar et al. (2005), is lower than unity.
Concerning the shear-induced contribution to the particle pressure, Boyer et al. (2011) proposed the following empirical formula for the dependance of sediment concentration α on the viscous number I v : www.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017 where B φ = 1 is a parameter of the dilatancy law (Boyer et al., 2011).Inverting Eq. ( 43) and substituting the definition of the inertial number I v gives the following expression for the shear-induced pressure: With Eq. ( 44), the total particle pressure p a can be calculated by Eq. ( 24).The frictional viscosity is defined by Eq. ( 41) with µ(I ) substituted by the µ(I v ) as shown in Eq. ( 42).

Grain inertia regime
In the grain inertia regime, the friction coefficient depends on the inertial number I = ∇u a d √ ρ a / p a and is calculated as with d the particle diameter, µ s the static friction coefficient, µ 2 an empirical dynamical coefficient, and I 0 an empirical constant of the rheology.For glass beads in air the typical values are as follows: µ s = 0.38, µ 2 = 0.64, and I 0 = 0.3 (Jop et al., 2006).
Concerning the shear-induced contribution to the particle pressure, it can be obtained from the dilatancy law α(I ) as proposed by Boyer et al. (2011) for the viscous regime of the granular-flow rheology.The adaptation to the inertial regime leads to the expression suggested by Maurin et al. (2016): According to Maurin et al. (2016), B φ = 0.31 in turbulent bed-load transport of spherical particles and Chauchat (2017) have proposed a value of B φ = 2/3 for sheet flows of nonspherical particles.
Inverting Eq. ( 46) and substituting the definition of the inertial number I gives the following expression for the shearinduced pressure (Chauchat, 2017): Similar to those in the viscous regime, the total particle pressure p a can be calculated by Eq. ( 24) and the frictional viscosity is defined by Eq. ( 41) with µ(I ) substituted by the µ(I ) as Eq. ( 45).
The rheology has been originally stated for steady uniform granular flows and this shear-induced pressure term induces a very strong coupling between the wall normal and the streamwise components of the particle-phase momentum balance equation.The granular-flow rheology has been used to simulate with success the transient flows such as the granular column collapse configuration (Lagrée et al., 2011) but the simulations have been performed at fixed volume concentration meaning that the shear-induced pressure term was neglected.In order to stabilize the model, a relaxation is added: old with a relaxation factor relaxPa that can be modified from the file constant/granularRheologyProperties.This is equivalent to assuming a relaxation in time for the dilatancy effect and the granular material does not dilate instantaneously to an imposed shear rate, which is physically justified.
The different closure laws implemented in SedFoam for the dense granular-flow rheology are summarized in Table 7.

Numerical implementation
The numerical implementation of the present Eulerian twophase flow sediment transport model is based on an opensource finite volume CFD library called OpenFOAM.Taking advantage of the numerical discretization schemes and framework of finite volume solvers in OpenFOAM, the twophase flow governing equations are implemented by modifying the solver twoPhaseEulerFoam (Rusche, 2002;Weller, 2002;Peltola, 2009).OpenFOAM uses the finite volume method over a collocated grid arrangement.The collocated arrangement stores all dependent variables at the cell center and the same control volume (CV) is used for all variables to minimize the computational effort.The advantage of the finite volume method is that the system of partial differential equations can be discretized on arbitrary three-dimensional structured or unstructured meshes.Thus, complex geometries can be easily handled.The Gauss theorem is applied to the convection and diffusion terms leading to conservative schemes.
To illustrate the numerical discretization, the fluid-phase momentum equation is taken as an example.Rearranging the fluid-phase momentum equation (Eq.4) by dividing βρ b , the resulting equation can be written as The last term in the above equation, the gradient of fluidphase shear stress, can be written according to Eq. ( 7) and expanded as follows: In the above equation, the first two terms on the RHS are treated implicitly while the last two terms are treated explicitly.By substituting the expanded shear stress formulation in the momentum equation, the following equation is obtained: It is more convenient to rewrite the above equation into a matrix form: The matrix A b is composed of the diagonal terms of the algebraic system associated with Eq. ( 50), whereas H b includes the off-diagonal terms and the source terms.R b is composed of the explicit drag term, the turbulent suspension term, the gravity term, and the explicit diffusion terms: The same process can be carried out for the solid-phase momentum Eq. ( 3) that leads to the following: where α at the denominator is substituted by α = α+α Small to avoid dividing by zero when the solid-phase volume concentration vanishes.This partial differential system of equations can also be written as a matrix equation as follows: where the source term R a contains the following terms: Following Rusche ( 2002) the terms involving the ratio of particle-phase volume gradient to the volume concentration are treated on the cell face level in the predictor-corrector algorithm.However, we noted the exception of the particlephase normal stress p ff gradient for which a reconstruction of the surface normal gradient at the cell center allows us to get more stable solutions.
The advantage of separating the RHS of the momentum equations as the sum of two terms, R and H, is for writing the velocity-pressure algorithm.A similar method as the Rhie and Chow (1983) one can be applied for the gradient terms.The details of the velocity-pressure algorithm are presented in the next section.

Velocity-pressure algorithm
The velocity-pressure coupling and the consequent oscillations in the pressure fields are resolved by using the Rhie and Chow method (Rhie and Chow, 1983).The PISO (Pressure Implicit with Splitting of Operator) algorithm is used to solve fluid and particle velocities (Rusche, 2002;Weller, 2002;Peltola, 2009).
First, the intermediate velocities (u a * , u b * ) are computed using the corresponding momentum equations (equations without the pressure gradient term): where A a −1 and A b −1 represent the inverse matrices of A a and A b , respectively.These intermediate velocities do not satisfy the mass conservation Eqs. ( 2) and (1).To enforce mass conservation for each phase, the pressure equation is constructed by considering the continuity equations for the mixture obtained as the summation of Eqs. ( 2) and (1): where the subscript f denotes variables interpolated at the cell faces.The two expressions shown above are equivalent by using the Gauss theorem.At the discrete level, this equation is written as where a f = u a | f .n|f S f and b f = u b | f .n|f S f denote the fluid-and particle-phase velocity fluxes at the cell faces, respectively, and S f is the cell face area associated with face f .In the present model, the method of Rhie and Chow ( 1983) is adopted in order to avoid velocity-pressure decoupling and oscillations.The velocity correction equations are written as at the cell faces, and In order to simplify the notations, the intermediate velocity and face fluxes are denoted as u a/b and a/b f , respectively, corresponding to the first two terms on the RHS of Eq. ( 59).The volume-averaged velocity, U * , and the corresponding averaged flux, * f , at the predictor step are defined as Taking the divergence of the volume-averaged mixture velocity given by the velocity correction Eq. ( 59) and imposing the incompressibility constraint, ∇ •U = ∇ •(αu a +βu b ) = 0, one can build the pressure equation as a function of the pre-dicted velocity or predicted face fluxes: The Poisson equation for the pressure is then written in discretized form at the cell face level as This equation leads to a matrix system written at the cell faces.The resulting algebraic system is usually solved using a multigrid solver (GAMG).The resulting pressure field p * is used for the correction step in which the fluid-and particlephase velocities and face fluxes are corrected using Eq. ( 59): The volume-averaged flux is also corrected according to * * In order to ensure the mass conservation an iterative procedure of N cycles is sometimes required.From our experience, three iterations (N = 3) is usually enough for a convergence.The finite volume discretization of the equations has not been shown here but all the details can be found in Jasak (1996) and Rusche (2002).

Summary of the solution procedure
The numerical solution procedure for the proposed twophase flow model is outlined as follows: 1. solve for sediment concentration α, i.e., Eq. ( 1 55) and ( 52); c. calculate u a * and u b * using Eq. ( 56) without fluid pressure gradient term; d. construct and solve the pressure Eq. ( 61); e. correct fluid and particle velocities after solving pressure and update fluxes Eqs. ( 62)-( 63); f. go to (a-e) if the number of loops is smaller than N (no tolerance criteria).
7. advance to the next time step.
In the above solution procedure, the velocity-pressure coupling steps are looped for N times.The advantage of this loop is to avoid velocity-pressure decoupling caused by the direct solving method.From our numerical practices, the loop number N = 1 to 3 is usually enough to give reasonably accurate results, and shows good convergence, especially for steady flows.The time step, t, can be set to a constant value or adjusted automatically based on two Courant numbers, one related to the local flow velocity and the local grid size Co = 1/2 f f t/V p (the same as for single-phase problems) and one related to the relative velocity Co r = 1/2 f a f − b f t/V p which is specific to the coupling of the fluid and sediment-phase momentum equations in the two-phase flow model.The most limiting time step is used as the criterion for setting the adjustable time step.Our practice is to set these two Courant numbers to 0.3 and 0.1, respectively.

Model verification and benchmarking
In this section, four benchmarking cases are presented to validate and verify the numerical implementation of the model.The first one concerns the pure sedimentation of a suspension of non-cohesive spherical particles for which experimental data are available.The second one concerns the laminar bedload problem for which an analytical solution exists.The third case is unidirectional turbulent sheet flow and the fourth one is about the scour at an apron.For each case, an analytical solution, experimental data, or empirical formula is used to validate the model.All the input files for the four cases described in the paper are available in the tutorial folder of the distribution and Python scripts based on an open-source postprocessing toolbox are available as well for facilitating the training on using SedFoam-2.0.

Pure sedimentation
The first test case corresponds to a pure sedimentation of non-cohesive particles; this test case allows us to validate the implementation of the pressure-velocity coupling algorithm when the flow is induced by the sediment phase.The other component that is tested here is the permanent contact pressure model (Eq.25) for p ff that allows us to predict a stable deposited sediment bed.The experimental dataset from Pham Van Bang et al. ( 2008 = 950 kg m −3 .The suspension is initially well-mixed in a cylindrical container (base diameter 50 mm, height 100 mm) with an initial solid volume concentration (α 0 = 0.5).The averaged concentration profiles are measured using a proton MRI device (Ecole de Ponts ParisTech, Champs-sur-Marne, France) with a vertical resolution of about 1 mm and a temporal frequency of 0.16 Hz (see Pham Van Bang et al., 2008, for details).This test case has been used by Chauchat et al. (2013) to validate a 1-D vertical two-phase flow model.A script for computing the solution is provided in the tutorial folder coming with the release 2.0 of SedFoam-2.0.
The mesh is composed of 200 cells in the vertical direction with a uniform distribution over a height h 0 = 0.06 m and the time step, t, is set to 10 −3 s.The numerical schemes for temporal and spatial derivatives are listed in Table 8, and details of these numerical schemes can be found in the user guide of OpenFOAM.The lateral boundaries are set to cyclic while the front and back boundaries are set to empty (i.e., 2-D problem).At the top boundary, the pressure is fixed with a zero value, and a zero gradient (i.e., Neumann boundary conditions) is imposed on all the other quantities.At the bottom, the velocity of both phases are set to zero while a fixedFlux-Pressure condition is imposed for the pressure.As mentioned by Lee et al. (2016), the fixedFluxPressure boundary condition must be used to avoid spurious oscillation on the wall normal pressure gradient close to the boundary.The granular rheology is turned on and an effective fluid viscosity model from Boyer et al. (2011) is used.
Figure 1 shows the comparison of the numerical results with the experimental data in terms of settling curves  = max{y | α ≥ 0.5 (α max +α 0 )}.The agreement between the numerical simulation results and the experiments is very good, suggesting that the closures for the drag and the particle pressure allows us to reproduce a pure sedimentation problem.This is the most basic test case for a two-phase flow sediment transport model.

Laminar bed load
The second test case is inspired by Chauchat and Médale (2010) in which the analytical solution for laminar bed load driven by a Poiseuille flow from Ouriemi et al. (2009) has been used to verify a three-dimensional numerical model.The goal of this test case is to verify the numerical implementation of the granular rheology in SedFoam-2.0against an analytical solution.The analytical solution can be derived only when using the Coulomb rheology (see Table 7) and the Einstein mixture fluid viscosity model (see Table 2).Those parameterizations will be used first for verification, then the numerical solution using the µ(I ) rheology will be compared with the previous numerical solution.The novelty compared with Chauchat and Médale (2010) is that the solid-phase concentration is obtained as a result of the model in SedFoam-2.0;however, it was imposed as constant in our earlier work.The details concerning the analytical solution can be found in Ouriemi et al. (2009) and will not be further detailed herein.A script for computing the solution is provided in the tutorial folder included with the release of SedFoam-2.0.
The numerical domain setup is based on the experimental configuration of Aussillous et al. (2013).The channel height is h 0 = 0.065 m and the particles are made of PMMA with a density ρ a = 1190 kg m −3 and a diameter d = 2 × 10 −3 m.The fluid density is ρ b = 1070 kg m −3 and the kinematic viscosity is ν b = 2.52 × 10 −4 m 2 s −1 .The pressure gradient is fixed to f x = 100 kg m −2 s −2 using the parameter gradP-MEAN in the file constant/forceProperties.The vertical domain is discretized into 200 uniform cells, and the time step is t = 10 −3 s.The numerical schemes are identical to the pure sedimentation case in Sect.4.1 (see Table 8).The lateral boundaries are set to cyclic while the front and back boundaries are set to empty (i.e., 2-D problem).The velocity of both phases is set to zero at the top and bottom boundaries while the pressure is fixed to zero at the top boundary and a fixedFluxPressure condition is imposed at the bottom boundary.Figure 2a, b, c show the comparison of the abovementioned analytical solution with the numerical solution in terms of sediment concentration (panels a, d), velocity (panels b, e), and particle pressure profiles (panels c, f).In the analytical solution, the sediment concentration profile is a step function with no particles in the upper half of the domain and with the maximum packing concentration in the lower half.The two-phase numerical model, based on continuous assumptions, is not able to reproduce exactly this sharp sediment concentration transition; this is because the sediment concentration profile is obtained using the momentum balance between gravity and the permanent contact contribution to the particle pressure (Eq.25).Despite this slight discrepancy, the numerical solution in terms of velocity profiles is in very good agreement with the analytical solution.Because the granular-phase viscosity is directly related to the particle pressure (Eq.40), the key issue for the granular rheology model is in the accurate prediction of the particle pressure profile.The comparison presented in the right panel shows that even if the agreement in sediment concentration profile is not perfect, the particle pressure profile is very close to the analytic solution.This explains the very good numerical prediction of the velocity profile.In Fig. 2d-f, the exact same problem is solved using the non-linear dense granular-flow rheology µ(I ) (see Table 7).This solution is compared with the solution of the 1-D code described in Aussillous et al. (2013).The very good agreement between the two numerical solutions allows us to validate the implementation of the µ(I ) rheology in SedFoam-2.0.It should be pointed out that, for this problem, the shear induced pressure contribution is not turned on, explaining why the sediment concentration and the particle pressure profiles are not affected by the change in the granular rheology.

Turbulent sheet flows
In this subsection the model results are compared with experimental results from Revil-Baudard et al. (2015) and Sumer et al. (1996) for turbulent sheet flows; the goal of these test cases is to validate the numerical implementation of different turbulence models (mixing length and k − ω) and to calibrate the free parameter B.
The first case is based on the experimental configuration of Revil-Baudard et al. (2015) for unidirectional sheet flows.The particles are non-spherical lightweight PMMA particles with density ρ a = 1190 kg m −3 and diameter d = 3 ± 0.5 mm.The measured settling velocity is W fall = 0.056 m s −1 .The fluid is water with density ρ b = 1000 kg m −3 and kinematic viscosity ν b = 10 −6 m 2 s −1 .The water depth is h f = 0.17 m, the energy slope is S f = 0.19 %, and the mean velocity is U = 0.52 m s −1 .
In the numerical configuration, the flow is driven by a pressure gradient (f x = 18.639 kg m −2 s −2 ).The mesh is composed of 400 cells in the vertical direction with a uniform distribution and the time step is set to t = 2×10 −4 s.Again, the numerical schemes are identical to the ones listed in Tawww.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017  9.
ble 8.The lateral boundaries are set to cyclic while the front and back boundaries are set to empty (i.e., 2-D problem).At the top boundary, the pressure boundary condition is fixed-Value or homogeneous Dirichlet and a zeroGradient or homogeneous Neumann is imposed on all the other quantities.At the bottom, the velocity of both phases are fixedValue or homogeneous Dirichlet, a fixedFluxPressure condition is imposed for the pressure and a zeroGradient or homogeneous Neumann boundary condition is used for the TKE and the TKE dissipation (ε and ω).The fixedFluxPressure is equivalent to set the second-order derivative of the pressure in the wall normal direction to zero.The results are presented in Fig. 3 for volume-averaged velocity profiles (panels a, d), defined as u = αu a + βu b , sediment concentration profile (panels b, e), and Reynolds shear stress profile (panels c, f).The numerical results are compared with the measurements reported in Revil-Baudard et al. (2015).Different combinations of granular stress model and turbulence model are presented.In panels (a, b, c) the µ(I ) rheology is used in combination with the mixing length (ML), k − ω, and k − ε models, while in panels (d, e, f) the kinetic theory model (KT) is used in conjunction with the k−ω and k−ε turbulence models.Note that the mixing length model has not been used with the kinetic theory as Eq. ( 36) requires an estimation of the TKE which is not straightforward to estimate when using a mixing length model.Amongst these different configurations, the µ(I ) rheology coupled with the mixing length model (panels a, b, c) corresponds to the model proposed by Revil-Baudard and Chauchat (2013) and Chauchat (2017), the kinetic theory coupled with the k − ε (panels d, e, f) corresponds to the model proposed by Hsu et al. (2004) and Cheng et al. (2017a).The mixing length model has been calibrated based on Revil-Baudard et al. (2015) data using the µ(I ) rheology (Chauchat, 2017) and the von Karman constant is reduced to κ = 0.225 to model a significant turbulence damping under sheet-flow conditions.Also, the rheological parameters of the µ(I ) rheology are identical to the ones proposed by Chauchat (2017) and they are summarized in Table 9.
In the k, ε, and ω equations, the turbulence damping is related to the drag term that is controlled by t mf , which parameterize the correlation between particles and fluid fluctuating motions t mf = e −B St (t mf = 1 at low Stokes numbers and t mf = 0 at high Stokes numbers).Based on Cheng et al. (2017b) large eddy simulations, the parameter B should be taken at 0.25 when using the kinetic theory.The same value Table 9. Physical parameters for the numerical simulations of Sumer et al. (1996) experiments.Amoudry (2014).After this calibration of the turbulence models, the results obtained with the different combinations of granular stress and turbulence models are discussed.In terms of velocity profiles, both the µ(I ) rheology and the kinetic theory are able to reproduce reasonably well the measurements provided that the turbulence model is adequately tuned.However, when using the k−ω or k−ε model with the µ(I ) rheology, the velocity profiles are underestimated suggesting that the dissipation in the sheet-flow layer is too strong and the velocity gradients are too small.The B value could be recalibrated to give a better result but this is not the purpose of the present contribution.A general conclusion is that the numerical solution does not depend too much on the choice of the turbulence model in between the k−ω and k−ε models.Concerning the sediment concentration profile, none of the two granular stress models are able to recover the measurements in the denser part of the sheet-flow layer.According to Revil-Baudard et al. (2015), this discrepancy might be related to the very strong near-bed intermittency and the assumptions used in the present turbulence-averaged formulation may be too simple.It may require a 3-D turbulence-resolving numerical simulation approach to capture this feature (Cheng et al., 2017b).

Case
In order to further assess the model, the same combinations of granular stress and turbulence models are applied to Sumer et al. (1996)  = 10 −6 m 2 s −1 .The physical parameter for the three cases SUM A, SUM B, and SUM C are summarized in Table 9.
The results in term of volume-averaged velocity profiles, defined as U = αu a + βu b , concentration, and shear stress profiles are presented in Fig. 4. Using the µ(I ) rheology with the mixing length model leads to an underestimation of the velocities and the velocity gradient in the sheet layer, at least for the first two cases.This could be explained by a smaller von Karman constant.The results obtained using the µ(I ) rheology with the k − ω and k − ε turbulence models again give very similar results that are in quite good agreement with the measured profiles.When using the kinetic theory and the k − ε model, the velocity profiles and the velocity gradients are overestimated in the sheet layer.Concerning the concentration profiles, they are all very similar and the results are not very sensitive to the different combinations of granular stress and turbulence models.For the shear stress profiles, it is observed that using the µ(I ) rheology, all the profiles are very close to one another.However, when using the kinetic theory, the Reynolds shear stress penetrates deeper into the sheet layer than when using the µ(I ) rheology.This is probably due to the presence of the fluid-particle interaction term in the granular temperature equation (Eq.36).
Different hypotheses can be proposed to explain the discrepancies presented above.According to Maurin et al. ( 2016), the µ(I ) rheology is able to describe accurately the dense granular-flow regime in turbulent sheet flows but it fails to predict the granular shear stress in the more dilute suspended layer at intermediate concentrations (α ≤ 0.3).It is therefore possible that the turbulence model is tuned in a way that the turbulent stress is overestimated leading to an underestimate of the streamwise velocity.Concerning the kinetic theory, it is well-known that it only applies when particle collisions are binary (Jenkins, 2006), i.e., when the concentration is not too high (α ≤ 0.3).For higher concentration, the particle-particle collisions involve more than two particles and the particle-particle contacts become "chattering".Jenkins (2006) proposed an extended kinetic theory that can tackle this problem.Moreover, the theory proposed by Berzi (2011) and Berzi and Fraccarollo (2013) has been applied with success to study turbulent sheet flows.Furthermore, we believe that the frictional stress model used with the kinetic theory is too simple and can not reproduce the granular behavior in the denser part of the sheet layer whereas the µ(I ) rheology is better suited.This is also supported by experimental observations from Capart and Fraccarollo (2011), the authors showed that the thickness of the dense frictional layer increases with the Shields number.A way to improve the kinetic theory would be to incorporate the latest theoretical developments proposed by Jenkins (2006) as done by Berzi (2011) and Berzi and Fraccarollo (2013) but this is beyond the scope of the present contribution.

Scour at an apron
In order to demonstrate the multidimensional capability of SedFoam-2.0, the fourth test case corresponding to the development of the scour downstream of an apron was examined.Following the numerical studies of Amoudry et al. (2008) and Cheng et al. (2017a), the problem is simplified as presented in Fig. 5 The top boundary is set as a symmetry plane.In all the simulations in this case, the numerical schemes are also similar to the ones in Table 8, except that "Gauss limitedLinear 1" is used for the divergence operators.As initial condition, the velocity of both phases, the sediment concentration, the TKE, and the TKE dissipation variables (ε or ω) are set based on one-dimensional simulation results using funkySetFields.The details of the boundary conditions are summarized in Table 10.The rough wall log-law is written as   11.
where u * = 3.69 cm s −1 is the bed friction velocity, κ = 0.41 is the von Karman constant, and k s = 2.5 d is the Nikuradse roughness length.Four combinations of fluid turbulence models (k − ε and k − ω) and granular stress models (KT and µ(I )) have been used (see Table 11 for reference).Figure 6 shows three snapshots of sediment concentration contours at different instants during the scour process, using k − ε and kinetic theory (left panels) and k − ω and µ(I ) granular rheology (right panels).At t = 10 s, the development of a scour hole near the inlet can be identified (see Fig. 6; top panels).As time elapses, the maximum scour depth increases and the scour perturbation is propagating downward.The snapshots presented in Fig. 6 show that results are qualitatively in agreement with one another using the different closures.According to experimental studies (e.g., Breusers, 1967;Breusers and Raudkivi, 1991), the development of the scour hole is rapid at the initial stage, and eventually reaches an equilibrium state.Breusers (1967) has suggested an empirical law to describe the rapid initial development of the scour hole: where T s is a characteristic timescale, and the exponent n s characterizes the speed of the scour development.Notice that Eq. ( 65) only describes the initial development of the scour depth, and the equilibrium scour depth can not be determined from this empirical formula.As the scour depth increases, the flow velocity reduces near the sediment bed, when the flow becomes weak enough; thus when it is below the criteria for sediment motion, an equilibrium scour depth can be obtained.The equilibrium scour shape is generally independent of the flow velocity and grain size if the Shields parameter is sufficiently large compared with the critical Shields parameter (Laursen, 1952;Chane, 1984).The development of the upstream bed angle, ζ s0 , can reach an equilibrium more rapidly.Breusers (1967) proposed the following empirical formula: where T ζ s0 is the equilibrium timescale for upstream bed angle and ζ ∞ s0 is the upstream bed angle at equilibrium.In Fig. 7, the numerical results for the four simulations presented above are shown in terms of these two quantities together with a best fit of the two empirical formula Eqs. ( 65) and (66).A summary of the fitted parameters is given in Table 11.First, the model is able to reproduce the power law for the initial development of the scour depth with values of n s in the range reported by other studies with higher T s values.The fitted values for the upstream bed angle are also in agreement with former studies.
In conclusion, this test case shows good capability of the proposed two-phase flow model to deal with multidimensional flow configurations.Further work is needed to improve the model validation and the model sensitivity to flow turbulence and rheological parameters.This requires more detailed experimental data that, to the best of our knowledge, are not available at present.

Conclusions
In this paper, a comprehensive two-phase flow model for sediment transport applications has been presented and the details concerning its implementation in OpenFOAM have been given.The proposed model provides different options for the modeling of flow turbulence (mixing length, k − ε, or k − ω) and intergranular stress (kinetic theory of granular flows or dense granular-flow rheology).The first validation test case presented on sedimentation of monodispersive spherical suspension allows us to validate the numerical implementation of the pressure velocity coupling and the modeling of the permanent contact contribution to the particulate pressure.The implementation of the dense granular-flow rheology, the mixing length and the two-phase k − ω models are original contributions.The dense granular-flow rheology is implemented using a regularization technique and is verified against an analytical solution for the laminar bed-load problem.The application of the model to turbulent sheet flows allows us to discuss the sensitivity of the model results to different combinations of intergranular stress and turbulence models.A first set of tuning coefficient values is provided, and the results are in reasonable agreement with four different experimental configurations.The last application on scour allows us to illustrate the multidimensional capabilities of the solver.The scaling laws proposed by earlier works are recovered by the model but the results are also sensitive to the modeling choices on the granular and turbulence models.In light of these model applications, some questions remain on the optimum values of the turbulence model coefficients, which will need more high-resolution measurements, for a wide range of flow conditions.The open-source numerical model presented here is expected to facilitate this endeavor in the future.
As a general conclusion, the aim of this contribution is to provide a comprehensive two-phase flow sediment transport modeling framework to the scientific community.Intense efforts have been made to ensure its reliability and numerical robustness.This numerical tool is suitable to address various physical problems for which the classical sediment transport approach is not working very well or requires more model assumptions.However, the readers are reminded that twophase flow simulations are still relatively time-consuming and require finer spatial resolution and smaller time steps than classical sediment transport models.
Appendix A: Notations

k
− ε model Cheng et al. (2017a) have implemented the k − ε model refined from Hsu et al. (2004) and Yu et al. (2010), in which www.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017 the volume concentration of fluid: β = 1 − α; 3. update the drag parameter K in the drag term, e.g., Eq. (5); 4. solve for the fluid turbulence closure, update k, ε, or ω (depends on the turbulence closure k − ε or k − ω), and then calculate the eddy viscosity and effective fluid total viscosity; Geosci.Model Dev., 10, 4367-4392, 2017 www.geosci-model-dev.net/10/4367/2017/ 5. solve for the particle-phase stress (kinetic theory model or the dense granular rheology); 6. PISO-loop, solving velocity-pressure coupling for N loops: a. construct the coefficient matrices A a and A b and explicit arrays H a and H b using Eqs.(54) and (51); b. update the other explicit source terms R a and R b , Eqs. ( ) is used for this validation.The suspension consists of mono-dispersed spherical polystyrene beads of diameter d = 0.29 ± 0.03 mm, and of density ρ a = 1050 kg m −3 in Rhodorsil silicone oil of viscosity ν b = 2.01 × 10 −5 m 2 s −1 and of density ρ b (panel a) and sediment concentration profiles (panel b).The settling curves show the time evolution of the vertical position of the upper Y up i and lower Y low i sediment interfaces.The upper one represents the transition between the suspension at the initial volume concentration and the clear water above it, while the lower one denotes the interface between the suspension and the granular bed at maximum packing fraction.They are defined as follows: Y up i = max{y | α ≥ www.geosci-model-dev.net/10/4367/2017/Geosci.Model Dev., 10, 4367-4392, 2017

Figure 1 .
Figure 1.Comparison of two-phase flow model results with experiments of Pham Van Bang et al. (2008) (a) settling curves: time evolution of the lower and upper interface positions (circles: experiments; lines: model) and (b) profiles of sediment concentration (dashed blue lines: experiment; solid red lines: model).

Figure 2 .
Figure2.(a-c) Comparison of the streamwise velocity profiles for the flow of a Newtonian fluid over a granular bed having a Coulomb rheology between two infinite parallel planes obtained by numerical simulations with the analytical solution ofOuriemi et al. (2009) in terms of sediment concentration (a, d), velocity profiles (b, e), and particle pressure (c, f) profiles.In panels (d-f) the same flow conditions are used but the granular rheology is the µ(I v) rheology and the results are compared with a one-dimensional numerical solution presented inAussillous et al. (2013).

Figure 3 .
Figure 3.Comparison of two-phase numerical results with experiments of Revil-Baudard et al. (2015) (black symbols) in terms of volumeaveraged velocity profiles (black lines) in the left column of panels, sediment concentration (red lines) in the center column of panels, and Reynolds shear stress (blue lines) and granular stress (red lines) in the right column of panels, using the dense granular-flow rheology (µ(I )) with the three turbulence models (mixing length, ML; k − ω; and k − ε) in (a, b, c) and the kinetic theory of granular flows (KT) with the three turbulence models (mixing length, ML, k − ω, and k − ε) in (d, e, f).The values of the different numerical and physical parameters are reported in Table9.

Figure 4 .
Figure 4. Comparison of two-phase numerical results with experiments ofSumer et al. (1996) (black symbols) in terms of volume-averaged velocity profiles (black lines) in the left column of panels, sediment concentration (red lines) in the center column of panels, and Reynolds shear stress (blue lines) and granular stress (red lines) in the right column of panels, using the dense granular-flow rheology (µ(I )) with the three turbulence models (mixing length, ML; k − ω; and k − ε) and the kinetic theory of granular flows with the k − ε turbulence model.The upper row of panels corresponds to the SUM A case, the middle row of panels corresponds to the SUM B case, and the bottom row of panels corresponds to the SUM C case.The values of the different numerical and physical parameters are reported in Table9.
experimental configurations corresponding to Shields numbers in the range θ ∈ [1.38; 2.1].For this configuration, the particles are made of acrylic, density ρ a = 1140 kg m −3 , and are of cylindrical shape with a mean diameter d = 2.6 mm.The measured settling velocity is W fall = 0.073 m s −1 .The fluid is water, density ρ b = 1000 kg m −3 and kinematic viscosity ν b

Figure 5 .
Figure 5. Sketch of the scour downstream of an apron.

Figure 6 .
Figure 6.Sediment concentration contour at different times during the scour process using k − ε and kinetic theory (left column of panels) and k − ω and µ(I ) granular rheology (right column of panels).
. The sediment bed is made of sand, density ρ a = 2650 kg m −3 and diameter d = 0.25 × 10 −3 m.The fluid is water with density ρ b = 1000 kg m −3 and kinematic viscosity ν b = 10 −6 m 2 s −1 .The flow depth is fixed to h 0 = 0.15 m, and the initial bed depth is h b = 0.05 m.The length of the domain downstream of the apron is L x = 1 m.A uniform grid is used in the streamwise direction ( x = 10 −3 m) while the mesh is refined vertically at the bed interface ( y ∈ [1.15×10 −4 −1.15×10 −2 ] m in the water column and y ∈ [1.15 × 10 −4 − 4.66 × 10 −4 ] m in the sediment bed).The bottom boundary, the lower part of the inlet (forming the step) and of the outlet are set as wall boundaries.The upper part of the inlet is an inlet boundary where the velocity profile is imposed according to the rough wall log law (Eq.64) and turbulent quantities are imposed as con-stant values following recommendation from the ESI group 1 .At the outlet, a directionMixed boundary condition is used for the velocities of both phases (zeroGradient or Neumann boundary condition for the streamwise component and fixed-Value or homogeneous Dirichlet boundary condition for the vertical component) and the hydrostatic pressure is imposed.

Figure 7 .
Figure 7. Numerical results for the temporal evolution of the normalized maximum scour depth and upstream bed angle.Different lines represent the best-fit curves for each run for which the parameters are given in Table11.
the rheology (s −1 ) D conductivity of granular temperature (kg m −1 s −1 ) f i external force that drives the flow (kg m −2 s −2 ) g s0 radial distribution function for dense rigid spherical particle gases (-) the granular rheology (-) J int energy dissipation (or production) due to the interaction with the carrier fluid phase (-)k turbulent kinetic energy (m 2 s −2 ) k s Nikuradse roughness length (m) K drag parameter (kg m −3 s −1 ) n s characteristic exponent (-) p fluid pressure (Pa) p aparticle normal stress (Pa) p a collisional component of the particle pressure (Pa) p ff permanent contact component of the particle pressure (Pa) q j flux of granular temperature (kg s −1 ) equilibrium timescale for upstream bed angle (s) t mf turbulent drag parameter (-) t l characteristic timescale of energetic eddies (s) t p particle response time (s) u k i phase k velocity, i = 1; 2; 3 represents streamwise, spanwise, and vertical components, respectively (m s −1 ) u * bed friction velocity (m s −1 ) U mean velocity (m s −1 ) W fall settling velocity (m s −1 ) Y low i position of the lower sediment interface for pure sedimentation test case (m) Y up i position of the upper sediment interface for pure sedimentation test case (m) α solid-phase sediment concentration (-) α 0 initial solid sediment concentration (-) α ∞ s0 upstream bed angle at equilibrium in the scour at an apron case (degrees) β fluid-phase volume concentration (-) δ s scour depth (m) t time step (s) ε turbulent kinetic energy dissipation rate (m 2 s −3 ) γ energy dissipation rate due to inelastic collision (kg m −1 s −3 ) κ von Karman constant (-) λ bulk viscosity (kg m −1 s −1 ) bed angle in the scour at an apron case (degrees)TableA2.List of the main model input parameters together with their default values when relevant.

Table 1 .
Summary of Eulerian two-phase models for sediment transport applications.

Table 2 .
List of the mixture viscosity model that can be selected through the keyword FluidViscosityModel together with model equation references.

Table 3 .
Default coefficient values for the k − ε model.

Table 4 .
Default coefficient values for k − ω model.

Table 5 .
List of the turbulence models that can be selected through the keyword RASModel.

Table 6 .
List of the kinetic theory closures that can be selected through the different keywords listed in the table together with model equation references.

Table 7 .
List of the dense granular-flow rheology models that can be selected through the keyword FrictionModel and PPressureModel together with model equation references.

Table 8 .
Numerical schemes for each term in the momentum and mass conservation equations used in the validation test Sect.4.1.ξ denotes a dummy variable for illustration purposes.

Table 10 .
Summary of the boundary conditions implemented in the 2-D scour downstream of an apron configuration.The following abbreviations are used: zG is zeroGradient for Neumann boundary condition, fV is fixedValue for Dirichlet boundary conditions, dM is directionMixed for mixed Dirichlet-Neumann boundary conditions, fFP is fixedFluxPressure corresponds to a Neumann boundary condition for the pressure gradient, and hp for Dirichlet boundary condition using the hydrostatic pressure.

Table 11 .
Summary of the numerical results obtained for the scour at an apron using the different combinations of turbulence and granular stress models and comparison with existing two-phase numerical results on this configuration.

Table A1 .
Table of notations.