Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using Smoothed Particle Hydrodynamics (SPH) accelerated with a Graphics Processing Unit (GPU)

Abstract A two-phase numerical model using Smoothed Particle Hydrodynamics (SPH) is applied to two-phase liquid-sediments flows. The absence of a mesh in SPH is ideal for interfacial and highly non-linear flows with changing fragmentation of the interface, mixing and resuspension. The rheology of sediment induced under rapid flows undergoes several states which are only partially described by previous research in SPH. This paper attempts to bridge the gap between the geotechnics, non-Newtonian and Newtonian flows by proposing a model that combines the yielding, shear and suspension layer which are needed to predict accurately the global erosion phenomena, from a hydrodynamics prospective. The numerical SPH scheme is based on the explicit treatment of both phases using Newtonian and the non-Newtonian Bingham-type Herschel-Bulkley-Papanastasiou constitutive model. This is supplemented by the Drucker-Prager yield criterion to predict the onset of yielding of the sediment surface and a concentration suspension model. The multi-phase model has been compared with experimental and 2-D reference numerical models for scour following a dry-bed dam break yielding satisfactory results and improvements over well-known SPH multi-phase models. With 3-D simulations requiring a large number of particles, the code is accelerated with a graphics processing unit (GPU) in the open-source DualSPHysics code. The implementation and optimisation of the code achieved a speed up of x58 over an optimised single thread serial code. A 3-D dam break over a non-cohesive erodible bed simulation with over 4 million particles yields close agreement with experimental scour and water surface profiles.


Introduction
Flows with two or more phases exhibit highly non-linear deformations and free surfaces are a common occurrence in applied hydrodynamic problems in mechanical, civil and nuclear engineering. The two-phase liquid-solid interaction is a typical problem in hydraulics and more specifically, flow-induced erosion. Other examples include port hydrodynamics and ship-induced scour, wave breaking in coastal applications and scour around structures in civil and environmental engineering flows.

A C C E P T E D M A N U S C R I P T
These subaqueous sediment scouring phenomena are induced by rapid flows creating shear forces at the surface of the sediment which causes the surface to yield and produce a shear layer of suspended particles at the interface and finally sediment suspension in the fluid. The current application is very difficult to treat with traditional mesh based Eulerian approaches due to the fluid-sediment interface, the highly non-linear deformation and fragmentation of the interface and the presence of a free surface leading to entrainment of the sediment particles by the liquid phase. These challenges require alternative simulation techniques. In the past two decades the novel Lagrangian approach Smoothed Particle Hydrodynamics (SPH) [12] has emerged as a meshless method ideal for this application. The scheme has been applied to a variety of problems such as free-surface flows [14], flood simulations [37], coastal flows [7], and geotechnical problems [3].
With numerous applications within engineering industries, there is a great deal of interest in non-Newtonian multi-phase flows. Rodriguez-Paz and Bonet [30] used the Generalised Viscoplastic Fluid (GVF) model to model the shear and plug flow of debris and avalanche failures as an non-Newtonian Bingham flow. Hosseini et al. [16] tested a variety of models such as the Bingham, power law and Herschel-Bulkley non-Newtonian models to examine non-Newtonian flows. Ran et al. [29] used a concentration based viscosity for the sediment phase similar to the work Shakibaeinia and Jin [32] for the Moving Particle Semi-Implicit (MPS) scheme. However, the aforementioned models examine mostly the rheological aspects of non-Newtonian flows. Recently Sibilla [33] applied the Exner equation to simulate the local scour caused by a 2-D horizontal wall jet on a non-cohesive granular bed downstream of a solid protection apron with reasonable success. Falappi et al. [8] adopted the Mohr-Coulomb criterion to model scour in reservoir flashing by using the Newtonian constitutive equation in a pseudo-Newtonian approach. Manenti et al. [20] compared the Mohr-Coulomb pseudo-Newtonian approach of Falappi et al. [8] with the Shield's criterion for the same experimental application. Ulrich et al. [36] used a similar approach and developed a scour and resuspension multi-phase model for ship-induced scouring near ports. Their model makes use of the Mohr-Coulomb approach for predicting the yielding of the sediment bed with a water-soil suspension model based on the Chezy-relation using piecewise linear interpolations between the soil, liquid and critical maximum viscosity for the suspension viscosity of the sediment. In a different approach, suited to geotechnics and embankment failures Bui et al. [2] replaced the simplistic plastic behaviour of the Mohr-Coulomb material by using an associated and non-associated plastic flow rule based on the Drucker-Prager model in combination with an elastic model of Hooke's law in the absence of the equation of state used previously by other researchers.
The aforementioned approaches tend to examine some flow features, in detail, but separately; i.e. the non-Newtonian characteristics of the sediment, the shear layer and the yielding independently. However, the rheology of sediment induced under rapid flows undergoes several states which are only partially described by previous research in SPH. This paper attempts to bridge the gap between the geotechnics, non-Newtonian and Newtonian flows encountered in scour by rapid liquid flows and sediment resuspension by proposing a model that combines the yield characteristics of sediment, the non-Newtonian rheology of the yielded sediment and a Newtonian formulation for the sediment entrained by the liquid applied to the SPH discretization scheme.
Note that our aim is to examine the rheology and scouring of the soil from a hydrodynamic approach and not from a geotechnical point of view. For more evolved geotechnics models the reader is referred to the work of Bui et al. [2].

A C C E P T E D M A N U S C R I P T
Adequately resolving the interface is essential to capturing complex industrial flows accurately with variable physical properties for each phase. Despite its suitability for such problems, SPH is well known for being an expensive method computationally [4]. Multiphase SPH simulations naturally involve many more particles further increasing the computational demands and costs [23]. In recent years, the massively parallel architecture of graphic processing units (GPUs) has emerged as viable approach to accelerate simulations requiring a large number of particles such as encountered in industrial applications. The Lagrangian nature of SPH makes the method not only ideal for large deformation flows with non-linear and fragmented interfacial multiple continua, but also makes it ideally suited to parallelisation on GPUs [4,15]. Accelerating SPH on a GPU is therefore the method of choice in this paper. Herein, we have modified the open source DualSPHysics solver [5] to include the two-phase liquid-solid model. DualSPHysics is a CPU/GPU weakly compressible solver package with pre-and post-processing tools capable of performing simulations on millions of particles using the GPU architecture targeted to real life engineering problems involving non-linear, fragmented and free-surface flows. This paper is organized as follows. Section 2 presents a brief description of the SPH discretization scheme. In Section 3, a detailed description of the numerical model including the sub-closure models and the multi-phase features such as the yield surface, constitutive modelling and sediment resuspension is presented. Speedup results from the GPU hardware acceleration are given in Section 4. Next, the numerical results and comparison with other numerical models and experimental results are shown for 2-D and 3-D cases, followed by the conclusions in Section 6.

SPH formalism
The basic principle of the SPH formulation is the integral representation of a function f which may represent a numerical or physical variable defined over a domain of interest Ω at a point x. The integral approximation or kernel approximation according to [12] reads with h defined as the smoothing length that characterizes the size of the support domain of the kernel and W the weighting or kernel function. The kernel function is chosen to be a smooth, isotropic and an even function with compact support (i.e. the finite radius of influence around x). In this paper, the fifth-order Wendland kernel with compact support of 2h has been used [40]: where the normalisation constant a d is 3/4h, 7/4hπ 2 and 21/16hπ 3 in 1-D, 2-D and 3-D space, respectively.
In a discrete domain Equation (2.1) can be approximated by using an SPH summation in the form of where V is the volume of the particle expressed as the ratio of the mass m to density ρ and N is the number of particles within the support. Throughout this paper the Latin subscript i denotes the interpolating particle and j refers to the neighbouring particles. The ... symbol denotes an SPH interpolation and will be dropped for simplicity in the rest of the paper. By dropping the approximation parentheses and the order of approximation term, the final form of the particle approximation in a discrete form is . More details of the SPH formulation can be found in [12] and more recently [40]. Figure 3.1 shows a schematic image of the physical domain and processes which are given in the following sections along with their SPH discretization. The domain consists of a saturated soil region which is subjected to the motion of fluid above which scours a region at the interface leading to a non-Newtonian flow in the yielded sediment. The interface between the yielded sediment and the un-yielded sediment is represented by the yield surface. At the interface, an approximation of Darcy's law forces accounts for fluid transfer across the yield surface only. Above the yielded region, sediment that has been suspended is assumed to be behaving according to a Newtonian fluid.

Navier Stokes equations and sub-closure models
The rheological characteristics of the domain are described by the Lagrangian form of the Navier Stokes equations discretized using the SPH scheme to approximate the multi-phase flows of this paper. Greek superscripts α, β denote coordinate directions by employing  (3.4) and

Equation of state
To approximate an incompressible fluid, the weakly compressible SPH approach (WCSPH) uses an equation of state (EOS) to link pressure to density in the form of (3.6) where ρ 0 is the reference density and and γ is the polytropic index with values between 1 to 7 and C s0 is the numerical speed of sound calculated as max 0 10u C s  . (3.8) where u max is the maximum velocity magnitude in the domain. Further information regarding the WCSPH approach can be found in [26].

Boundary conditions
The wall boundary condition applied in this paper is the dynamic boundary conditions (DBCs) [6] where particles representing the wall are organised in a staggered arrangement and satisfy the same equations as the fluid particles but their position and velocity are prescribed. The advantages of the DBC include the straightforward computational implementation and the treatment of arbitrary complex geometries. This makes them particularly amenable to use within the GPU code DualSPHysics. Boundary conditions discretization is not the focus of the present research; for more information the reader is referred to more recent work [10].

Time integration
The time stepping algorithm is an explicit second-order predictor-corrector integrator scheme. The scheme predicts the evolution in time at half time steps. These values are then corrected using the forces at half time steps, followed by the evaluation of the values at the end time step [13]. The scheme is bounded by the CFL condition, the maximum force term and the numerical speed of sound as demonstrated by Monaghan [24]. In addition, an extra restriction is imposed based on the viscous forces. The CFL condition reads, (3.9) where f i is the force per unit mass of particle i and Co is the Courant number set to 0.3 in this paper and ν is the kinematic viscosity. The same time integration scheme is used for both phases using the minimum Δt resulting from the CFL condition.

Newtonian viscous formulation
For a more consistent formulation of the multiple phases we write Stokes' theorem for a general fluid using the thermodynamic pressure and the extra stress tensor in the form of (3.10) where δ αβ is the Kronecker delta. The theorem assumes that the difference between the stress in a deforming fluid and the static equilibrium stress is given by the function f determined by the rate of deformation D. When f is linear for an isotropic material, such as water, the fluid is called Newtonian and the constitutive equation can be written in the general form as (3.11) where   denotes the strain rate tensor. For incompressible flow   = D αβ since the D γγ is zero by the continuity equation. In the WCSPH formalism Equation (3.11) can be used to obtain the total stresses in the momentum equation, thus the strain rate tensor of Equation (3.11) can be calculated for the velocity gradients as (3.12) Therefore, viscous stress tensor can be calculated from the Newtonian constitutive equation that relates the strain rates to the viscous stresses by (3.13) The total viscosity μ of Equation (3.13) represents the dynamic and eddy viscosity μ τ obtained through the Smagorinsky algebraic eddy viscosity model by (3.14) Herein, for the 3-D simulations a Large Eddy Simulation (LES) model has been used to model the turbulent characteristics of a multi-phase flow. The LES model is a standard Smagorinsky algebraic eddy viscosity model as described by Dalrymple and Rogers [7] in a WCSPH formalism for Newtonian fluids. The model reduces to a mixing length model in 2-D simulations and at low velocity flow regimes reduces to a laminar viscosity model. This study has not considered if the turbulence has been resolved sufficiently, however it is an important aspect of the liquid-sediment flow field and should be investigated further. On a different approach, such as of the Shields criterion the soil particles are resuspended mostly from the interface [20]. The turbulence characteristics are of primary importance for the suspension of the sediment soil within such a model by considering mean local variables.

δ-SPH
In this paper the δ-SPH approach is used in the fluid phase and sediment phase independently. Thus, the computation of the δ-SPH term for the fluid does not include sediment and vice versa. δ-SPH accounts for the bulk viscosity dissipation in the mean pressure by a dissipation term applied to the continuity equation similar to the Stokes condition. However, the current δ-SPH formulation is based on an empirical artificial dissipation which is not related to the bulk viscosity, in a similar manner to the artificial viscosity of Monaghan [27]. The continuity equation in WCSPH does not guarantee a divergence-free velocity field and the equation takes the form of (3.17) where δ d is a parameter usually set to 0.1 [21]. The first term on the right-hand side of Equation (3.15) is the velocity divergence and the second term is the δ-SPH dissipation to the continuity equation. A comprehensive description of the implementation can be found in [5].

Particle shifting
While particles' field variables such as velocity and pressure are well predicted with WCSPH, issues may arise in negligible dynamics or large dynamics respectively [25].
Both states occur often in multi-phase flows and are a frequently addressed issue. Manenti et al. [20] and recently Ulrich et al. [36] used the XPSH approach of Monaghan [25] with an additional smoothing applied to the particle position through a smoothed velocity. Additionally, Ulrich et al. [36] used the XSPH formulation to perform pressure smoothing for large dynamics flows where the pressure exceeded twice the hydrostatic pressure. However, such smoothing procedures tend to smooth the dynamics of the system such as sharp interfaces and discontinuities.
In this work, using the particle shifting algorithm of Lind et al. [18] liquid or yielded sediment particles are shifted a very small distance each time step from areas of high particle concentration to low concentration by using a Fickian-type approach to maintain a more regular particle arrangement (3.18) where C is the particle concentration and n and s is the vector normal and tangent unit vector respectively. The D' parameter is the diffusion coefficient of Skillen et al. [34] based on a von Neumann stability analysis using a constraint based on the velocity magnitude of the particle that reads , (3.19) Finally, b is the bi-tangent unit vector to account for shifting at 3-D free surfaces and interfacial flows that has been recently applied by Mokos et al. [22] to WCSPH.
The surface treatment extension to 3-D of Mokos et al. [22] is only applied to the liquid phase since most large dynamics are dominant in the liquid phase.

Sediment model
The saturated sediment rheological characteristics induced by the liquid flow field exhibit different behavioural regimes shown in Figure 3.1 that adhere to the sediment properties and shear stress of the liquid phase at the interface. The non-Newtonian nature of sediment flows results from several physical processes including the Mohr-Coulomb shear stress τ mc , the cohesive yield strength τ c which accounts for the cohesive nature of fine sediment, the viscous shear stress τ v which accounts for the fluid particle viscosity, the turbulent shear stress of the sediment particle τ t and the dispersive stress τ d which accounts for the collision of larger fraction granulates. The total shear stress can be expressed as (3.20)  Yielding: At low stress state the sediment remains un-yielded in that region with the yield strength of the material being greater than the induced stress by the liquid phase and is dominated by the first two terms on right-hand side of Equation (3.20). Nevertheless, the saturated sediment stress state should be taken into accountedsee Section 3.4.1.  Constitutive modelling: In a high stress state the sediment is yielded and behaves as a non-Newtonian rate dependant fluid using the last three terms on right-hand side of Equation (3.20). Herein, the dispersive stresses are not modelled explicitly as with other discrete sediment models. Instead the constitutive formulation accounts for these dispersive stresses [30]. Typically sediment behaves as a shear thinning material with a low and high shear stress state of a pseudo-Newtonian and plastic viscous regime respectively [17] see Section 3.4.3.  In addition, an approximation of the generalised Darcy law has been applied in order to simulate the saturated soil motion and the interaction of the sediment and water at the interface and the saturated yielded sediment phase described in Section 3.4.5.

Yield surface
To determine the state of the sediment SPH particle (yielded or un-yielded region), a yield criterion is used to relate the maximum shear strength of the soil sediment to the hydrodynamic shear strain at the fluid-soil interface. Above a predefined value related to the shear strength, the sediment is assumed to be at rest whereas below the critical threshold the sediment undergoes yielding. Note that compression is assumed to be positive throughout this article.
Considering a simple shear case where no motion in the sediment phase takes place until a critical value of shear stress τ y is reached, the fluid stresses acting on the sediment are in equilibrium with the yield strength of the sediment [9] i.e.
where τ y is defined as the sum of the Coulomb and cohesive yield strength as (3.22) and J 2 is the second invariant of the deviatoric shear stress tensor τ αβ defined as (3.23) Recalling Equation (3.11), the rate dependant isotropic Newtonian fluid expression for the viscous stresses is written as (3.24) Squaring both sides of Equation (3.11) the following equality is derived (3. 25) where the term II D is the second invariant of the strain rate tensor defined as (3. 26) Thus, the critical threshold for the sediment yielding at the interface can be written as ( 3.27) At this point, a yield criterion for the sediment phase is needed to provide the critical value of the sediment shear stress. In this study, the Drucker-Prager yield criterion has been used following previous investigation by Fourtakas et al. [9] indicating the suitability of different yield criteria.
The Drucker-Prager model is written in a general form as The parameters a and κ can be determined by projecting the Drucker-Prager (DP) onto the Mohr-Coulomb yield criterion in a deviatoric plane. The yield surface of the deviatoric plane at π/6 corresponds to where φ is the internal angle of friction and c the cohesion of the material.
Finally, using Equation (3.27) yielding will occur when the following equation is satisfied Although, the Drucker-Prager criterion or other Coulomb based yield criteria are simplistic and depend on the pressure, internal angle of friction and cohesion of the sediment, they have been used extensively in SPH ([1], [9], [20] and [36]) with satisfactory results.

Sediment skeleton and pore-water pressure
Equation (3.27) assumes a constant critical value of shear stress τ y . This might not always be true for saturated drained conditions. Sediment pressure changes according to the lithostatic conditions and the pore water pressure for a fully saturated sediment. In isotropic, fully saturated sediment under drained conditions the Terzaghi relationship holds where subscripts t, eff and pw denote the total, effective (skeleton) and pore-water pressure, respectively.
The total pressure of the system can be calculated simply by accounting for the hydrostatic and saturated pressures as (3.32) where h is the height, γ is the unit weight and subscripts w and sat denote the water and saturated phase respectively as shown in the schematic of  Equation (3.32) requires the surface to be tracked in order to determine the maximum height of the saturated sediment which is usually computationally expensive [20]. Instead, the equation of state can be used by modifying the reference pressure dependant on the numerical speed of sound of Equation (3.7) by relating the pore water pressure to the saturated sediment pressure as (3.33) where B is based on the fluid properties where the subscript pw refers to the pore water, sat to the saturated sediment and w to the water phase, thus recovering the pore water pressure in the saturated sediment even though the density ratio is still based on the saturated sediment.
The total pressure in the sediment is calculated by using Equation (3.6). The sediment skeleton (or effective) pressure can finally be calculated using Equation (3.31). Note that the skeleton pressure can only be applied to fully saturated soils. A partly saturated sediment methodology can be found in [36].

Constitutive modelling
The rheology of the shear mobile layer of the sediment at the interface can be described using viscoplastic rheological laws usually described by Bingham models [17,30]. The Bingham model is one of the simplest models and provides a satisfactory description of the viscoplastic behaviour of subaqueous sediment flows but it cannot approximate the range of stress regimes encountered in scouring, i.e. pre-and post-yield behaviour. Nevertheless, a variety of other Bingham models such as the bi-viscosity and Herschel-Bulkley (HB) models are often used in subaqueous flows mimicking the Bingham rheology of a viscoplastic material in low and high stress states [17]. However, using the HB model poses numerical issues for zero shear stress states. To avoid this issue, the Herschel-Bulkley-Papanastasiou (HBP) model [28] has been employed to model the rheological characteristics of the yielded region. The HBP model reads (3.35) where m controls the exponential growth of stress, n is the power law index and μ is the apparent dynamic viscosity. Figure 3.3(a) shows the initial rapid growth of stress by varying m whereas Figure 3.3(b) shows the effect of the power law index n.
The advantage of the HBP model is the pseudo-Newtonian region defined by the growth of stress parameter m and the power law index n in the plastic region. This two region approach in combination with the yield criterion has been chosen to model the soil phase without the use of an explicit elastic branch. Nevertheless, successful elastoplastic models have been developed in SPH and applied to sediment transport [2]. The HBP model provides information on the pre-yielded and post-yield regions after the apparent yield region defined by the Drucker-Prager criterion with a low stress and high stress region. The sediment phase can also be modelled as a typical shear thinning fine grained material. In addition, there is no need for scale-back methods as used in previous work by other researchers [36].
For a specific skeleton pressure, the inequality of Equation (3.27) defines the yielded surface at that point (or particle). Regardless of whether the particle is yielded or not the shear stress is calculated using the HBP model with the specific yield stress. However, in the un-yielded region the sediment particles motion is restricted by setting du/dt = 0 but discontinuities in the stress summations of the momentum Equation (3.17) do not arise since the viscous stresses in the un-yielded region are computed and assigned to the sediment particle. For the suspended entrained sediment particles a concentration suspension viscosity (see Section 3.4.5) is used to avoid "particle freezing" and force imbalances [36].

Approximation of Seepage forces in the yield surface
To simulate the saturated soil motion, the interaction of the sediment and water phases should be taken into account. The behaviour of saturated soils is determined by the interaction between the soil skeleton and the pore water pressure. When the mixture is deformed the sediment skeleton is compressed and pore-liquid flows though the pores. Water seeping through the pores of a soil produces a drag force on the sediment phase originating from viscous forces. This force acts on the direction of the water flow [3]. Darcy's law is often used to impose the viscous drag force as (3.36) where subscripts w and s denote water and sediment, respectively, K is based on the soil characteristics and can be written as (3.37) where n r is the porosity and k is the soil permeability and γ = ρg is the unit weight. The seepage can be added in the momentum equation as an extra term and Equation (3.2) becomes For simplicity, it is assumed that the water does not flow in the un-yielded region and seepage only acts at the interface of the un-yieldedyielded regions and the interface (see Figure 3.1). Also, the soil mixture is assumed to be isotropic and fully saturated under drained conditions. Although the assumption that water does not flow in the un-yielded region is not strictly correct for the accurate representation of the seepage forces in the soil body, Darcy law forces are approximated in the on the yield surface which is our point of interest with this article. Other SPH practitioners ( [3], [31]) modelled seepage forces using Darcy law, by using a two SPH particles layers approach to superposition the liquid and soil layer. Unfortunately, this technique tends to be cumbersome and memory intensive. GPUs are memory restricted and such a 3-D model would not be feasible with the current technology.
In the SPH formalism the seepage force can be expressed in SPH form as where the term 0.01h 2 in the denominator is included to avoid singularity as x ij → 0. The seepage force is only applied to particle i in the yielded region for all j particles irrespective of the phase. A similar methodology has been proposed by [3].

Suspension
At the interface, the fluid flow at a sufficient large velocity (τ y << τ fluid ) will suspend the sediment particles in the fluid. This sediment entrainment by the fluid can be controlled by using a concentration volume fraction of the mixture in the form of where the summation is defined within the support of the kernel and j sat refers to the yielded saturated sediment particles. The size of the concentration sampling is chosen as to adhere with the kernel support size of SPH. When a sediment particle is suspended, it is modelled as a pseudo-Newtonian fluid using Equation assuming an isotropic material with spherically shaped sediment particles. Equation (3.42) is applied only when the volumetric concentration of the saturated sediment particle within the SPH kernel is lower than 0.3 which is the upper validity limit of Equation (3.42). Hence, when a yielded sediment particle volumetric concentration is below the threshold of 0.3 which coincides with the validity of the Vand equation, the sediment particle is treated as a Newtonian fluid, retaining its properties with the exception of the viscosity that follows the Vand Equation (3.42).

Hardware acceleration using GPUs
Recently, the introduction of GPU cards to scientific computing provided acceleration to massively data-parallel computations including n-body simulations. The attraction comes from the parallel, multithreaded many-core computing architecture of GPU cards which is well suited to problems of data-parallel computing since the same instance is executed on many threads in parallel. Such data-parallel computations in SPH are the particle interactions of the interpolated particles. NVIDIA GPU cards use an instruction set architecture, namely CUDA in C/C++ programming language to use the massive parallelism of GPUs. DualSPHysics, an SPH C++/CUDA solver, exploits the massive parallelism of GPU cards resulting in speedups comparable to large CPU clusters [4],and has been extended in this work for multi-phase water-sediment flows.
The multi-phase model described earlier has been implemented in the DualSPHysics CUDA code [5] in order to accelerate the computations and allow higher particle resolution. The higher resolution aims to capture the important yielding and viscous effects at the interface of the multi-phase model where variable viscosity and large density variations affect the shear layer of the suspended particles. In addition, extension to three dimensional (3-D) and real life industrial applications requires a large number of particles (on the order of millions). An in-depth investigation of different GPU techniques on implementing multi-phase models to GPU cards can be found in the recent work of Mokos et al. [23].
Mokos et al. [23] gas-liquid modelling involved only some extra terms in the governing equations and therefore minimum impact on memory. Nonetheless, as demonstrated by the authors branching is a major bottleneck in the GPU architecture, creating different lists of particles within CUDA kernels calls can partially resolve the branching impact on the code acceleration. However, memory occupancy and register usage remains unresolved.
The sediment phase rheological model differs from the Newtonian formulation considerably and includes extra calculation such as the yielding of the surface, the pore water pressure estimation, the seepage forces and the suspension of the sediment in the fluid under low concentration. These extra operations have to be accommodated in the CUDA kernels code and consequently, branching and memory occupation becomes a major concern in multiphase modelling due to the limited flow control and memory cache of the of the GPU architecture which is primarily based on arithmetic operations.
A different approach has been applied. Instead of calling the CUDA kernels depending on the phase of the particle, a memory operation is performed that points to the constants of the specific phase. These constants, such as the reference density, viscosity, polytropic index, yield strength, etc, are stored in the constant memory of the GPU, a fast access memory

A C C E P T E D M A N U S C R I P T
which is common to the streaming multiprocessors. This avoids branching and extra register usage but increases arithmetic operations. Most importantly, it does not impact the memory occupancy of the local memory and register space which is the foremost requirement of the multi-phase model. Nevertheless, since most operations are carried out despite the phase of the particle, and factoring the extra calculations of the sediment phase the arithmetic operations increase which is a different requirement to Mokos et al. [23] which was mostly bounded by branching. Figure 4 shows the speedup achieved in relation to a single core single thread C/C++ code running on an Intel i7 processor 2.6 GHz with the same multi-phase implementation in 2-D.
A satisfactory speedup of x58 is achieved for 1.5 million particles on a NVIDIA Tesla K20 GPU card.

Numerical results
The DualSPHysics code has been previously validated for single-phase non-linear flows [5]. This section presents numerical results, first for a validation case for yielding of the sediment phase separately, and then fully saturated conditions with both phases combined. Finally, results are presented for a simulation of a 3-D dam break over an erodible bed comparing with experimental data.

Droplet impact of a flat plate and the effect of particle shifting
To demonstrate the effectiveness of the particle shifting algorithm and the δ-SPH formulation under potentially violent impact flows a test case with a droplet impacting a horizontal plate is employed [19]. The radius of the 2-D sphere is 8.5 x 10 -4 m impacting on a flat plate with a particle spacing of 2 x 10 -5 m resulting in a total of 9835 particles. The liquid droplet has a reference density ρ l = 1680 kg/m 3 with a viscosity of μ l = 6.4 x 10 -3 Pa.s without gravity.
Results obtained by just using δ-SPH are compared with a simulation that uses δ-SPH and particle shifting. Figure 5.1 shows a comparison of the two configurations 370 s after impact and the effects of the particle shifting algorithm on the formation of unphysical voids.

A C C E P T E D M A N U S C R I P T
(a) δ-SPH only (b) particle shifting and δ-SPH The void structure of Figure 5.1 (b) is created due to particles following the streamlines after a sudden impact of the droplet to the plate creating particle line structures (stacks of particle lines following the streamline) that eventually collapse after particle clumping has occurred. This unphysical void and particle clumping phenomena disturbs the pressure field either by creating spurious pressures or by pressure oscillations manifested as voids at later time as shown in Figure 5.1 (b). In addition to the void formation effect of shifting, the pressure field of Figure 5.1 (a) is considerably smoother and symmetric. However, since particles are shifted from their original position, at impactkeeping in mind that in WCSPH small compressions of the order of 1-2% are permitted through the equation of state and the numerical speed of soundthe compressed particles are shifted from the high concentration impact zone to the lower concentration interior droplet domain.
This redistribution of particles from the high to the low concentration zone maintains the pressure field in the droplet transmitting the pressure wave and therefore numerically reducing the compressibility of the fluid. Figure 5.2b shows that using the particle shifting algorithm of Section 3.3.3 generates a pressure field inside the droplet that is in closer agreement with a Volume-of-Fluid (VoF) reference solution [19].  Also, the yield surface shows satisfactory agreement with the experimental data and the numerical results of Bui et al. [2]. Similar issues near the wall have been identified by other researchers using the dynamic boundary conditions of DualSPHysics in multi-phase simulations [22]. The repose angle of the current SPH model was found to be 16˚ which is 1˚ degree more that the more

A C C E P T E D M A N U S C R I P T
involved model by Bui et al. [2]. These results are satisfactory considering the simplistic nature of our soil model.

2-D Erodible dam break
Having shown example validation cases for the sediment phase, the paper now presents validation cases for two phases combined.
The shear layer and suspension of the sediment is qualitatively validated using an experimental erodible dam break [11]. A definition sketch is shown in Figure 5.4. When the water column is released, the hydrodynamic stress at the interface induces yielding and scouring at the sediment interface. The shear layer induced propagates to a depth where the viscosity, density and pressure change from their initialized value.
In the computational setup a fully saturated sediment bed with height H s = 0.6 m is located below a dam of liquid with height H l = 0.1 m and length L l = 2.0 m. The particle spacing is set to dx = 0.002 m producing 328 353 particles in the domain.   [36] and current model at t = 1.00 s.

A C C E P T E D M
A N U S C R I P T Figure 5.5 shows the liquid and sediment profile at t = 0.25 s which follows the initial dam release. The profile of the SPH simulation is in good agreement with the experiment data following a similar trend for the liquid surface and the liquid-soil interface. A root mean square (RMS) error has been used to quantify the deviation of the SPH interface location from the experimental profile results by using the absolute position of the interface profile. The RMS error for the liquid-soil interface deviated by 1.15%. A departure is visible at the toe front of the dam were the numerical model run off distance is marginally forward of the experimental. This can been seen clearly in Figure 5.9 where a zoom of the dam toe is provided. However, the liquid-sediment SPH model over predicts scouring at t = 0.75 s and onwards in Figure 5.7 and Figure 5.8. Although the SPH model prediction is wider, the maximum depth is predicted correctly with an RMS error of 0.89% for the SPH model. An explanation of the over prediction could be due to the rather simplistic modelling of turbulence and yield criterion. More elaborate models with a non-associated plastic flow rule such as Bui et al. [2] would be advantageous.
Concluding, the new model predicts the liquid surface and interface profile satisfactory. The scouring profile was predicted adequately for the dam toe at the beginning of the simulation at t = 0.25 and slightly over estimated after t = 0.75 s in the current model. The results are reasonable, given our assumptions of isotropic stress, turbulence and the modelling of the yield surface with the Drucker-Prager criterion.

A C C E P T E D M
A N U S C R I P T

3-D Erodible dam break
This section presents numerical results for a fully 3-D validation case for dam-break flows over mobile beds involving fast transient flows and sediment transport. The experimental test case of Soares-Frazão et al. [35] provides validation data for numerical models. In the numerical model the initial particle spacing was set to 0.005 m for both liquid and sediment resulting in 4 million particles. This is the finest resolution that can be simulated by using an NVIDIA K20 GPU, restricted by the memory size of the GPU. It resulted in 17 particles over the initial depth of sediment. The initial density of the fluid and sediment was set to hydrostatic and lithostatic conditions, respectively. The δ-SPH parameter for this simulation was set to 0.1 as recommended in [21]. The fluid dynamic viscosity was 0.001 Pa.s and sediment viscosity was set to 150 Pa.s with the HBP m and n parameter set to 100 and 1.8 respectively. The value for the exponential growth m parameter value was chosen to approximate a Bingham model as closely as possible with a minimal pseudo-Newtonian region and the power-law exponent n to resemble shear thinning materials as shown in Figure  3.3.
Finally a small amount of cohesion was given to the sediment phase of c = 100 Pa to stabilise the interface and control the scouring near the dam gate. Next, the water level height and sediment profiles are presented against the experimental data for the aforementioned control points.
A C C E P T E D M A N U S C R I P T

Sediment bed profiles
The sediment profile in comparison with the experimental data is shown in Figure 5.11 for three different cross-sections of the sediment bed at locations y 1 = 0.2 m, y 2 = 0.7 m and y 3 = 1.45 m. The numerical results are compared with four different experimental runs (labelled as b 1 , b 2 , b 3 and b 4 ) as reported by Soares-Frazão et al. [35] with the SPH data superimposed on the experimental data. The data from runs b 1b 4 show the variability in the experimental data.
The RMS error was calculated similarly with the 2-D case of Section 5.3. However, the average experimental profile was computed first as reference data. The bed profile at y 1 shows satisfactory agreement with the experimental results for most of the length of the bed with an RMS error of 1.57% compared with the average experimental profile. Some deviations from the experiment are notable specifically near the dam break gate around x = 0.5 m to x = 1.5 m. Nevertheless, there is some variation in the repeatability of the experimental results of an RMS error of ±1.33% with some of the runs having lower scouring at the front with a peak for runs b 2 and b 4 whereas the numerical results are in better agreement with runs b 1 and b 2 . Also a small deviation is observed at x = 4.0 m where the numerical model over predicts the scouring region.
At y 2 , the agreement is marginally deteriorates from y 1 with small deviations near the wall at x = 0.5 m which were expected due to the boundary conditions implemented in DualSPHysics that can exhibit a sticking behaviour near the walls and an RMS error of 1.82% of the average experimental profile. The RMS error of the experimental profile is ±1.57%.
Finally, the y 3 bed profile shows similar behaviour near the wall over predicting the sediment height. Also, the sediment peak is slightly under predicted with a small delay on the location of the peak. The model RMS error for the y 3 profile is 2.94% of the average experimental profile.
A C C E P T E D M A N U S C R I P T However, due to the complexity of the 3-D dam break case the sediment profiles are satisfactory in view of the challenge posed by the hydrodynamics of the fluid especially at the gate with a rarefaction wave and an initial hydraulic jump and the SPH wall boundary conditions.

Water level measurements
The hydrodynamics of the flow are linked to the sediment scour mechanisms. In this section, two water level probes locations are used to measure the numerical water levels and compare with the experimental results of the 3-D dam break. The experimental profiles for gauge US1 and US6 shown in Figure 5.10 are compared with the SPH results. The water height levels of gauge US1 which is located near the gate, shows reasonable agreement with the reported experimental results with a RMS error of 2.94% over the average experimental profile and a RMS error of ±1.34% for the three experimental runs. It is notable that the larger over prediction is at the end of the simulation t = 20 s.
Similar results are shown for gauge US6 with an RMS error of 3.23%. However the variation between the experimental runs is smaller with an RMS error of ±0.68% and therefore resulting in a less accurate prediction by the SPH model. Also notable is a small dip at t = 8 s. Both graphs show slightly over predicted water heights particularly at location US6. This corresponds with the sediment profiles which show a small bed thickness over prediction downstream.
Concluding, the results from the numerical experiment of a 3-D dam break have been compared with a benchmark dam break case of Soares-Frazão et al. [35]. To the authors best knowledge this is the first time this test case has been performed with SPH mainly due to the large domain and therefore the high computational cost. The 3-D erodible dam break took 14 A C C E P T E D M A N U S C R I P T days on a Tesla K20 GPU card for 4 million particles. The simulations would benefit from recent advances in variable resolution [38] to make the simulation more efficient.
The 3-D dam break bed profiles of the sediment located at the bottom of the tank were satisfactorily reproduced in this numerical experiment with about 2-3% deviations near the gate and a small over prediction downstream that might be accounted for by the departure in the dynamics of the fluid. Also the water level profiles at two discrete locations have been presented with reasonable agreement to the experimental results.

Conclusions
In this paper the development and validation of a Smoothed Particle Hydrodynamics (SPH) multi-phase model has been presented. The multi-phase model focuses on liquid-sediment flows and more specifically the scouring and resuspension of the solid phase by liquid induced rapid flows. The choice of modelling technique in this work is based on explicit treatment of the liquid and solid phase using a Newtonian and a non-Newtonian constitutive model respectively that is supplemented by the Drucker-Prager yield criterion to predict the yielding characteristics of the sediment surface and a suspension model based on the volumetric concentration of the sediment.
Graphic processing units (GPUs), with massively parallel capabilities have been the choice of hardware acceleration in this work. GPUs' parallel architecture is well suited to n-body simulations. The open-source weakly compressible SPH solver DualSPHysics was chosen as the platform for the GPU implementation and optimisation of the multi-phase model with a reported speedup of 58 that allowed the simulation of large domains with millions of particles such as the erodible dam break 3-D case which was unfeasible in the CPU implementation.
Comparison between experimental and 2-D numerical results showed reasonable agreement on the interfacial and liquid free-surface profiles. More specifically the yielded region of the 2-D erodible dam break was captured satisfactorily against the experimental results. Also the free-surface evolution of the dam was predicted adequately. In comparison with previous numerical models the results were satisfactory with some improvements on the yielding characteristics of the sediment. However some sediment transport over prediction occurred at later times of the simulation. A 3-D case was used to validate the 3-D numerical model that was accelerated using a GPU coprocessor. The numerical results were sufficiently close to the experimental data for the scouring profile. Small deviations from the experimental results were present, however, as seen from the repeatability of the experiments over different runs large deviations in the experimental data were also present. Finally, the free-surface and water-level elevation was well predicted with similar deviations as observed for the sediment scouring profile across the bed. However, care should be taken when choosing the model properties as it can lead to over prediction of the scouring and water level elevation profile.
A C C E P T E D M A N U S C R I P T

A C C E P T E D M A N U S C R I P T
Graphical abstract