## Abstract

Sediment scouring is a common example of highly dynamic sediment transport. Considering its complexities, the accurate prediction of such a highly dynamic multiphase granular flow system is a challenge for the traditional numerical techniques that rely on a mesh system. The mesh-free particle methods are a newer generation of numerical techniques with an inherent ability to deal with the deformations and fragmentations of a multiphase continuum. This study aims at developing and evaluating a multiphase mesh-free particle model based on the weakly compressible moving particle semi-implicit (WC-MPS) formulation for simulation of sediment scouring. The sediment material is considered as a non-Newtonian viscoplastic fluid, whose behavior is predicted using a regularized *μ(I)* rheological model in combination with pressure-dependent yield criteria. The model is first validated for a benchmark problem of viscoplastic Poiseuille flow. It is then applied and evaluated for the study of two classical sediment scouring cases. The results show that the high-velocity flow currents and the circulations can create a low-viscosity region on the surface of the sediment continuum. Comparing the numerical results with the experimental measurements shows a good accuracy in prediction of the sediment profile, especially the shape and dimensions of the scour hole.

## INTRODUCTION

Highly dynamic (or rapid) movement of sediment material is a problem of common occurrence in many hydro-environmental and engineering settings. Reservoir flushing, dam break flow over sediment bed, sediment scouring around structures (e.g., under pipes or around bridge piers), riverbank failure, and underwater landslides, are a few examples of the highly dynamic sediment movement, characterized by rapid motion of dense sediment mass with large deformations and fragmentation and frequent topology change. Among these examples, scouring has a special importance due to its effect on the stability of the hydraulic structures. The highly erosive flow (with high velocity and strong flow circulation) is often the driver of sediment scouring. The foundations of contemporary research on sediment processes (e.g., erosion, transport, deposition), can be found in many experimental, field, and soft-computing studies. Sediment scouring has also been well studied using these techniques (e.g., Azamathulla & Ghani 2010; Dogan & Arisoy 2015; Guan *et al.* 2016; Park *et al.* 2016; Roushangar *et al.* 2017; Barati *et al.* 2018).

With the advancement in computing powers and the numerical algorithms, the numerical simulations have gained much attention. Traditional approaches for simulation of sediment processes have mostly been based on single-phase modeling, in which the flow equations are coupled with the advection-diffusion equation (for suspended sediments) and the empirical or semi-empirical relationships (van Rijn 1984) for determining processes such as bed-load transport. These methods, which are the base of the majority of the existing models, have proved their effectiveness in some past studies (e.g., Falconer & Owens 1990; Kashyap *et al.* 2017).

Nevertheless, they are not suitable for dealing with the highly dynamic movement of sediments (the case of sediment scouring) as they are not able to predict the underlying physics in detail. Multiphase models, with separate water and sediment phases, have proven to be an enabler for many sediment transport simulations. These methods are based either on continuum-based mesh-based Eulerian or Lagrangian methods such as finite element method, finite volume method, and material point method (Hsu *et al.* 2003; Larese *et al.* 2013; Lee *et al.* 2016), or discrete-based Lagrangian methods such as discrete element method (Calantoni *et al.* 2004; Tao & Li 2018). The first group is ill suited in dealing with the multiphase flows with large interfacial deformations and fragmentations due to mesh-adaptability and connectivity problems (Shakibaeinia & Jin 2011) and the second group is computationally expensive for large-scale problems as they simulate each individual grain of sediment (Rycroft *et al.* 2010). The mesh-free particle (Lagrangian) methods developed for continuum mechanics problems, such as smoothed particle hydrodynamics (SPH) (Gingold & Monaghan 1977) and moving particle semi-implicit method (MPS) (Koshizuka & Oka 1996), provide simple and natural yet powerful tools for the simulation of highly dynamic multiphase flows. The mesh-free particle methods can accurately capture any interfacial deformations and fragmentations. Therefore, they can be ideal candidates for dealing with the complexities involved in highly dynamic sediment movements. They have the Lagrangian nature of discrete-based methods, yet do not have their scalability issue, as they model the multiphase continuum.

Particularly, the MPS, on which this study is based, has proven its capability for many complex fluid flow problems (e.g., Gotoh & Sakai 2006; Ataie-Ashtiani & Farhadi 2006; Shakibaeinia & Jin 2012a, 2012b; Gambaruto 2015). Shakibaeinia & Jin (2010) proposed a weakly compressible version of the MPS method (WC-MPS) to improve efficiency and also the numerical fluctuations associated with the original MPS. Jafari-Nodoushan *et al.* (2016) improved the accuracy of the model near the solid and inflow boundaries. Shakibaeinia & Jin (2012a) proposed a multiphase MPS method for multi-density multi-viscosity systems. Fadafan & Kermani (2017) improved pressures' stability properties with the MPS method. MPS is also very versatile in the simulation of non-Newtonian flows, by adapting various rheological relations. By combining the multiphase MPS model with visco-plastic rheological relations such as Herschel–Bulkley and Bingham plastic, Shakibaeinia & Jin (2011, 2012b) developed a model for simulation of flow-driven sediment flows for the example application of mobile-bed dam break and sand jets and plumes. Fu & Jin (2015) and Nabian & Farhadi (2016) used a similar approach for simulation gravity-driven sediment flow for the case of submerged landslides. Similarly, the SPH method was used for simulation gravity-driven (Capone *et al.* 2010) and flow-driven sediment flow (Manenti *et al.* 2011; Khanpour *et al.* 2016). More recently, Jafari-Nodoushan *et al.* (2018) and Tajnesaie *et al.* (2018) have developed and fully characterized a multiphase WC-MPS model for gravity-driven sediment flows. Unlike past studies, they were based on a pressure-dependent viscoplastic model with a dynamic inter-grain mechanical pressure.

Another important aspect of the continuum-based numerical models is the constitutive (rheological) model used to mathematically describe the non-Newtonian behavior of a granular (sediment) continuum. While the general visco-plastic models such as Bingham plastic, and Herschel–Bulkley have been two of the widely used models in past studies, more recent attention has been directed towards *μ*(*I*) rheological model (Jop *et al.* 2006). Unlike the Bingham plastic and Herschel–Bulkley models, the *μ*(*I*) rheological model is specialized for granular flows. In this model, the rheological properties are explicitly related to the physical properties of granular material (e.g., grain size and density), and the granular flow regimes (quasi-static, inertial, and viscous regimes) are taken into account. Furthermore, Tajnesaie *et al.* (2018) showed that the *μ*(*I*) model produces slightly more accurate results than the Herschel–Bulkley model.

The present research therefore aims at developing and implementing a mesh-free particle (Lagrangian) model based on multiphase WC-MPS formulation for the simulation of two-phase highly dynamic movement of sediment for the specific case of sediment scouring. The method treats the multiphase system as a single multi-viscosity and multi-density continuum. An exponentially regularized *μ*(*I*) rheological model with the fully dynamic inter-grain mechanical pressure (unlike past studies) is applied to predict the non-Newtonian behavior of the sediment phase. The benchmark case of viscoplastic Poiseuille flow is used to evaluate the capabilities of the model dealing with non-Newtonian viscoplastic behaviors. The model is then validated and evaluated for simulation of the sediment scour around a submerged pipe and scour resulting from water jets (in a plunge pool) with different angles.

## GOVERNING EQUATIONS

*t*is time,

**r = r**(

*x, y*) is the position vector,

**u = u**(

*u, v*) is velocity vector,

*ρ*is density,

**g**is gravity acceleration,

*p*is pressure, and

**is deviatoric stress (shear stress) tensor given by: where**

*τ***I**is the unit tensor,

*η*is the effective viscosity of materials (given by the rheological model),

*ξ*is the second coefficient of viscosity, and

**E**is the rate of strain tensor. Here, in this study, the flow is considered weakly compressible and an equation of state (EOS) predicts the pressure using the density field

*p*

*=*

*f*(

*ρ*). For a weakly compressible fluid, the divergence of the velocity vector goes to zero , therefore the stress tensor reduces to

**2**

*τ*=*η*

**E**. Considering the two-phase system of sediments and water as a single multi-density multi-viscosity continuum, one set of governing equations (Equation (1)) is used for both phases (similar to Shakibaeinia & Jin (2012a) and Hu & Adams (2007)) and no additional interfacial forces will be required.

### Sediment rheology

*η*. While the effective viscosity of water phase is constant (independent of time and stress) that of sediment phase is determined through a rheological model (constitutive relationship). Among the rheological models and theories developed for granular and sediment flows, a frictional model based on the

*μ*(

*I*) law (Jop

*et al.*2006; Forterre & Pouliquen 2008) has recently gained wide attention. The

*μ(I)*rheological model is derived from the analysis of a large number of laboratory and numerical data, providing a general framework for describing a wide range of granular material behavior. Before defining the

*μ*(

*I*) model, one should remember the required basic definitions of tensor calculus, in which for

*D*-dimensional arbitrary tensor

**A =**[

*a*], the trace is given by , transposed by

_{kl}**A**

*, the first invariant by*

^{T}*A*, second invariant by

_{I}*A*, and magnitude by . The

_{II}*μ*(

*I*) model is a viscoplastic model, in which the material behaves as a rigid body for stresses less than a yield stress,

*τ*and flows as a viscous fluid for stresses larger than

_{y}*τ*. This model defines the effective viscosity as: where

_{y}*μ(I)*is the inertial-dependent coefficient of friction, and

*p*

_{s}is mechanical pressure (normal stress) between the sediment grains. Here, the yield stress,

*τ*is calculated from the Drucker–Prager yield criterion (Chen & Ling 1996), which for cohesionless material is calculated as

_{y}*τ*=

_{y}*μ*

_{1}

*p*where

_{s}*μ*

_{1}represents the first friction coefficient (a tangent of angle of repose

*θ*). The coefficient of friction

_{r}*μ*(

*I*) depends on a dimensionless number, called the inertia number

*I*, defined as the ratio between microscopic and macroscopic time scales, and depending on the regimes of granular flow. The coefficient of friction is given by (Jop

*et al.*2006): where

*μ*

_{2}is the second friction coefficient (

*μ*

_{1}and

*μ*

_{2}represent the upper and lower limits of

*μ*(

*I*)),

*I*is a reference inertial number (a constant coefficient, commonly determined through experiments). For a submerged granular flow with a low-viscosity ambient fluid (Stokes number =

_{0}*d*(

_{s}*ρ*)

_{s}p_{s}^{0.5}/

*η*>> 1) the inertial number is given by (Forterre & Pouliquen 2008): where

_{f}*d*is sediment grain size and

_{s}*c*is the drag coefficient, and

_{d}*ρ*and

_{f}*η*are the density and viscosity of the ambient fluid, respectively. Therefore, the effective viscosity can be calculated by:

_{f}*et al.*2018). Nevertheless, it relates the post-failure rheology to a variable term, which not only is a function of the material properties, but also depends on the pressure and strain rate. Equation (6) is singular and discontinues when the strain rate goes to zero (

**E**→ 0), therefore it cannot be directly used in the numerical model. To avoid this issue, a regularized continuous form has to be used. Here, we combine it with a widely used exponential regularization proposed by Papanastasiou (1987) to Equation (7) as: where

*m*is an exponent that controls the viscosity growth. This equation is valid for stresses smaller and greater than the yields stress and is continuous as for

**E**→ 0 the effective viscosity goes to the finite maximum value of

**.**

## NUMERICAL METHOD

*r*. A weight function (kernel function),

_{e}*W*(

*r*,

_{ij}*r*), determines the effect of a neighboring particle

_{e}*j*(with position vector

**r**

*) on the target particle*

_{j}*i*(with position vector

**r**

*), where*

_{i}*r*= |

_{ij}**r**

_{j}**−r**

*| is distance between two particles. This study uses the third order polynomial kernel function, proposed by Shakibaeinia & Jin (2010) as: A dimensionless parameter called particle number density,*

_{i}*n*, is used as a measure of density as: here 〈 − 〉 is the weight averaging operator. The MPS approximation of the spatial derivatives (interpolation, gradient, divergence, and Laplacian) for the arbitrary scaler

*f*and vector

**f**on target particle

*i*is given by (Shakibaeinia & Jin 2010): where n

_{0}is the initial particle number density and

*λ*is a coefficient defined by:

### Pressure gradient

*c*is an artificial (numerical) sound speed, selected to keep the compressibility (density change) lower than an acceptable value (typically less than 1%). The pressure gradient is then approximated using a modified form of MPS gradient formula as: where

_{0}### Rheological model implementation

As was mentioned, the model of this study considered the two-phase sediment water systems as a single multi-density multi-viscosity continuum. The fluid-type and sediment-type particles are used to discretize the continuum, where the fluid-type particles represent ambient water, and the sediment-type particles represent sediment grains and pore water between the grains (Figure 2). Each particle is assigned a density and viscosity. The water density is assigned to the fluid-type particles, and sediment bulk density, *ρ _{b}*, is assigned to the sediment-type particles.

*Ω*

_{s}) particles, and the total particle number density as . This volume fraction is then used for determining the interface and calculation of the density field. The effective viscosity,

*η*, of each particle is also needed, in order to calculate the divergence of shear stress tensor (;

**2**

*τ*=*η*

**E)**. The effective viscosity of sediment-type particles is variable (a function of stress and strain rate) and obtained from the rheological model, while that of fluid-type particle is constant (

*η*

*=*

*η*). From the

_{f}*μ*(

*I*) rheology formulation, discussed before, one can calculate the effective viscosity of a sediment-type particle

*i*as: The MPS approximation of two-dimensional strain-rate tensor,

**E**, on particle

*i*is given by:

*p*, for incompressible dry granular flow is equal to the pressure

_{s}*p*, given by the EOS. However, for submerged sediments (the case of this study), the pore pressure,

*p*, must be subtracted. MPS method is known to be associated with some unphysical pressure fluctuations, which can affect yield behavior. To avoid this issue, some of the past studies (e.g., Shakibaeinia & Jin 2011; Manenti

_{pore}*et al.*2011) have used a static pressure (assuming the vertical acceleration to be small negligible), as: where

*h*

_{i}is vertical distance of particle

*i*and the particle at sediment–water interface. However, this study uses a more generalized dynamic pressure, proposed by Jafari-Nodoushan

*et al.*(2018), as: This formula accounts for the vertical accelerations, and minimizes the effect of pressure fluctuations as it uses the smoothed value of thermodynamic pressure, 〈

*p*〉

*.*

_{i}*i*and

*j*, with viscosities of

*η*and

_{i}*η*, are interacting, none of the viscosities can be used. Shakibaeinia & Jin (2011) proposed a harmonic mean of viscosities for such interaction, given by

_{j}*η*= 2

_{ij}*η*/(

_{j}η_{j}*η*+

_{i}*η*). Therefore, the MPS approximation of the divergence of shear stress tensor for particle

_{j}*i*can be written as:

### Boundary conditions

Free surface pressure is assigned to free surface boundaries. To locate the free surface particles, the value of particle number density is used. Particle number decreases while approaching the free surface. This is used as the criterion for recognition of free surface particles. A particle will be determined as a free-surface particle if its particle number density is less than *β* times of the initial average particle number density . In this paper, a *β* value of 0.94 (recommended by Shakibaeinia & Jin (2010)) is used.

For rigid boundaries (bed and walls), a few layers of so-called ghost particles are considered outside the boundary to assign the boundary values and to avoid unwanted density deficiency near the boundary. The pressure of ghost particles is extrapolated from the interior of the computational domain (to assure the boundary repulsive force). For the particles within the boundary layer, a logarithmic velocity profile is used. The details about the boundary treatment can be found in Jafari-Nodoushan *et al.* (2016).

### Solution method

*α*is chosen between 0 and 1 (here

*α*

*=*0.5). In each time step the rheological model is applied to calculate the effective viscosity, which is then used to predict the velocity. Once the predicted velocity is calculated, particles are displaced, then the predicted value of the particle number density

*n**is calculated. The pressure is then computed from this particle number density. The pressure gradient term is then recalculated and used for calculation of corrected velocity and particle position. To account for the effect of high Reynolds flow, especially in water phase, a turbulence model is applied to calculate an eddy viscosity which is then added to the effective viscosity. Here, a large eddy simulation sub-particle scale (SPS) model (Shao & Gotoh 2005) is used. Details about the MPS implication of this model are available (Shakibaeinia & Jin 2010; Jafari-Nodoushan

*et al.*2016). This model has been used because it has been well tested for particle methods. Time step size must satisfy the Courant–Friedrichs–Lewy (CFL) stability condition as well as stability of the diffusion equation as:

where *d _{p}* is distance between the particles (size of the particles) and

*C*is the Courant number and

*η*is maximum viscosity.

_{max}## TEST CASES

In this section, first the model is validated for a benchmark, a viscoplastic Poiseuille flow, to evaluate the accuracy of the model in simulation of yielded and non-yielded regions in the flow field. Then, the developed model is validated and evaluated for simulation of sediment scouring in two common scouring situations: (1) sediment scouring under pipes and (2) sediment scouring induced by a water jet.

### Viscoplastic Poiseuille flow

The geometry of the viscoplastic flow between two static parallel plates is shown in Figure 3(a). The flow is driven by the gravity, **g**. Note that flow occurs only if *ρ*|**g**| > *hτ _{y}*, where

*τ*is a constant yield stress. Here, a dimensionless case with parameter values of

_{y}*h*

*=*0.5,

**g =**(0, −4/3),

*τ*=

_{y}*μ*

_{1}×

*ps*

*=*1, has been considered. The rheological model is assumed to be linear (

*I*< <

*I*

_{0}) to simplify the analytical solution as of a Bingham-plastic fluid and the material properties have been set to have equal to 0.5. The regularization exponent is set to

*m*

*=*200. Domain is represented by particles with size of 0.0125. The steady-state analytical solution is taken from Chatzimina

*et al.*(2005). In order to visualize the deformations, the particles are initially colored to form several parallel horizontal layers. Figure 3(b) shows a snapshot of the numerical results for the particle position at

*t*

*=*2 (when time simulation reaches the steady condition). The particle layers show two distinguishable regions, yielded (near walls) and non-yielded (near center). The velocity numerical and analytical velocity profiles have also been compared showing a very good compatibility, with RSME of 0.008.

### Sediment scour under circular pipe

In this section, the local scour beneath a submerged circular pipeline, based on the experiment of Mao (1986), is modeled. The sediment materials have a grain size of , bulk density of , and internal friction angle of degrees. The deposited layer has a thickness of 0.25 m, the water depth was 0.35 m and the pipeline diameter was 0.1 m. A two-dimensional computational domain with a length of 45 D and a height of 6 D is used in the simulation. The distance between the pipeline center and the inlet boundary is 15 D. The depth-averaged velocity at the inlet boundary is specified as 0.65 m/s. Figure 4 shows the general schematic of the forms of motion for the sand grain beneath the pipe.

The computational domain is represented with 28,508 particles (with initial distance of 0.01 m). The computational time on a machine with an Intel Core™ i5 (2.9 GHz) processor and 8GB of memory is about 5 hours for 1 second of simulations. Figure 5 shows results of numerical snapshots of local scour around a circular pipe problem at *t**=* 1.5, 2.5, and 5 seconds. The high velocity flow under the pipe and flow circulation create a scour hole that grows as the time goes on, until an equilibrium is reached. The eroded materials are deposited further downstream and create a sediment ridge. Comparisons of the experimental and numerical snapshots show similar scouring and deposition features. A quantitative comparison of experimental and numerical sediment surface profiles (Figure 6) shows a good agreement, especially for the scouring zone. The dimensionless depth of scour hole (*D _{s}* =

*d*/

_{s}*D*) is considered as the measure of the scouring.

The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole. The numerical deposition zone, however, seems to be more shifted to upstream. This can be due to numerical assumption used for calculation of pore and mechanical pressure (which assumes a negligible exchange of water mass between ambient water and pore water). Figure 7 presents the time evolution of the dimensionless depth of scour hole. The experimental and numerical results are quite compatible. Figure 8 provides a sample snapshot of simulated, volume fraction, velocity magnitude, and velocity vectors at *t**=* 1.5 seconds. In this figure, the flow pattern, characterized by high-velocity flow under the pipe and circulation behind the pipe, can be seen. This erodes the sediment materials under the pipe and deposits them further downstream. The viscosity field differentiates the immobile particles (with very large viscosity) and particles in incipient motion (low-viscosity particles). These low-viscosity particles can be observed under the pipe, where the erosion happens. The present model has not used any empirical relation for calculation of erosion and deposition of sediments.

### Jet-induced sediment scouring

Outflow from most hydraulic structures is in the form of a jet, which can be a source of severe sediment scouring at downstream channels. Here, the sediment scouring problem due to a non-submerged water jet (with various jet impact angles) in a rectangular plunge pool is investigated with the proposed mesh-free model of this study. The experimental data, required for comparison purposes, are taken from Pagliara *et al.* (2008). Figure 9 shows schematic geometry of the problem and Table 1 shows the geometrical parameters for two setups (A and B).

Setup | h/_{0}D | Jet diameter D (m) | Angle α (°) | Discharge Q (L/s) |
---|---|---|---|---|

A | 5.1 | 0.03069 | 60 | 1.5 |

B | 7.2 | 0.03069 | 45 | 1.5 |

Setup | h/_{0}D | Jet diameter D (m) | Angle α (°) | Discharge Q (L/s) |
---|---|---|---|---|

A | 5.1 | 0.03069 | 60 | 1.5 |

B | 7.2 | 0.03069 | 45 | 1.5 |

The water jet in these selected setups for modeling purposes has impact angles of 45 and 60° (with respect to a horizontal reference line). The sediment materials have a grain size of , bulk density of , and internal friction angle of . The depth (*h _{0}*) and width of the channel and thickness of the sediment layer are 0.7 m, 0.5 m, and 0.3 m, respectively. Parameter is dimensionless time. Numerical modeling was performed for two setups, A and B, with jet impact angles of 45 and 60°. The solution domain was represented with particles with average distance (particle size) of 0.00115 m. Figure 10 compares the experimental and numerical snapshots (showing particle types) of jet scouring problem for 60° jet (setup A) at

*t*

*=*3, 7, 15, and 47 seconds. The jet flow and the resulting circulations create a high-shear stress region (whose shear stress exceeds the yield stress), where the sediment material starts to be eroded. As a result, a scour hole forms and grows as the time goes on, until equilibrium is reached. The eroded materials are deposited further downstream under the gravity force and create a sediment ridge. Similarly, Figure 11 shows numerical snapshots of 45° jet (setup B) at

*t*

*=*25, 50, 75, and 100 seconds. As expected, the 45° impact angle results in a less deep and more stretched scour hole compared to the 60° impact angle. Comparisons of the experimental and numerical snapshots show similar scouring and deposition features.

Figure 12 provides a quantitative comparison of experimental and numerical sediment surface profiles for the 60° jet (Case A). Note that for Case B, the experimental snapshots are not clear enough for accurate extraction of the profiles. Figure 12 shows a good compatibility of the profiles, especially for the scour hole. The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole.

Two parameters, including the dimensionless depth of scour hole (*D _{s}* =

*d*/

_{s}*D*) and dimensionless height of deposition ridge (

*D*=

_{d}*d*/

_{d}*D*), are considered as the measures of the scouring and deposition features. Figure 13 presents the time evolution of

*D*and

_{s}*D*for jet angles (

_{d}*α*

*=*60 and 45°). The experimental and numerical results are quite compatible (especially for the scour hole depth) for both jet angles. The numerical model slightly underestimates the deposition ridge height as it stretches the deposition region to the downstream. The root mean square error provided in Figure 14 also shows a better prediction of the scour hole depth compared to the deposition ridge height. Accuracy of the results for the sediment deposition can be due to the numerical assumption used for calculation of pore and mechanical pressure (which is not suitable for suspended materials), as well the two-dimensional simulation of a three-dimensional case.

Figure 15 provides a sample snapshot of simulated viscosity, velocity magnitude/vector, and volume fraction at *t**=* 5 s. The flow circulation pattern has also been highlighted in this figure. As the jet flow hits the sediment surface, it is divided into two branches: (1) a spoon-like counterclockwise strong flow at the direction of the jet and (2) a weak backward (clockwise) flow in the opposite direction of the jet. These flow currents erode the sediment material, creating a scour hole, and deposit them downstream and upstream, respectively. As time goes on, the scour hole grows. Where the shear stress due to these flow currents exceeds the yield stress (i.e., low viscosity regions) the sediment materials start to move and separate from the bed. The material is carried by the flow current and deposited further downstream under their weight. This creates an active (yielded) sediment layer on the sediment surface. The sediment layer under this active layer remains relatively stagnant (un-yielded). This is because of the relatively low shear stress, high yield stress, and effective viscosity (due to the weight of the above layers).

## CONCLUSIONS

Multiphase mesh-free particle modeling of local sediment scouring with *μ*(*I*) rheology was developed for continuum modeling and evaluated for simulation of sediment scouring. The sediment material is considered as a non-Newtonian viscoplastic fluid whose behavior is predicted using a regularized *μ(I)* rheological model in combination with pressure-dependent yield criteria. The model is first validated for a benchmark problem of viscoplastic Poiseuille flow. Then, the developed model is validated and evaluated for two classical sediment scouring cases, namely, scouring under submerged pipes, and scouring resulting from water jets with different angles. The proposed numerical model has been able to successfully reproduce the shape and size of the scour hole. The numerical deposition zone, however, seems to be more shifted to upstream. This can be due to the numerical assumption used for the calculation of pore and mechanical pressure. The experimental and numerical results for the jet are quite compatible (especially, for the scour hole depth) for both jet angles. The numerical model slightly underestimates the deposition ridge height as it stretches the deposition region to the downstream. The results show how the highly erosive flow and circulations can create a low-viscosity region at the surface of the sediment continuum, where the shear stress can exceed the yield stress and the sediment can be eroded to form the scour hole. The results of this study not only demonstrated the capability of mesh-free particle methods for sediment scouring modeling, but also provided a better understanding of flow and sediment scouring in the test cases of this study.