Two-Phase Flow Modeling for Bed Erosion by a Plane Jet Impingement

: This paper presents experimental and numerical studies on the erosion of a horizontal granular bed by a two-dimensional plane vertical impinging jet to predict the eroded craters’ size scaling (depth and width). The simulations help understand the microscopic processes that govern erosion in this complex ﬂow. A modiﬁed jet-bed distance, accounting for the plane jet virtual origin, is successfully used to obtain a unique relationship between the crater size and a local Shields parameter. This work develops a two-phase ﬂow numerical model to reproduce the experimental results. The numerical techniques are based on a ﬁnite volume formulation to approximate spatial derivatives, a projection technique to calculate the pressure and velocity for each phase, and a staggered grid to avoid spurious oscillations. Different options for the sediment’s solid-to-liquid transition during erosion are proposed, tested, and discussed. One model is based on uniﬁed equations of continuum mechanics, others on modiﬁed closure equations for viscosity or momentum transfer. A good agreement between the numerical solutions and the experimental measurements is obtained.


Introduction
Vertical jet-induced scour erosion of soil has been studied for many industrial applications, such as aerospace or hydraulic engineering [1][2][3][4][5][6]. As reported by Metzger et al. [1,2], in aerospace engineering, soil erosion and crater formation could generate problems for launching and landing spacecraft. In hydraulic engineering, we have several examples, such as the work of Rouse [3] for testing criteria on erosion, the study of Hanson and Cook [4] assessing in situ the erodibility of soil material, and the research developed by Hanson and Hunt [5] in the so-called jet erosion test (JET). We also have the work of Perng and Capart [6] for water injection dredging (WID) and trenching in harbors by jet translation and possible jet inclination.
Based on experiments of 3D axisymmetric jets impinging on a smooth non-erodible wall, Rajaratnam and Beltaos [8] and Phares et al. [18] divided the flow into three zones: • The free jet region. There is no influence of the boundary. It can be described by self-similar solutions [19,20]. • The impinging region. It is characterized by a stagnation point on the jet axis at the wall with a deflexion of flow streamlines [21]. • The wall jet region. It is close to the wall but away from the impinging zone where the flow is mainly parallel to the wall [22][23][24].
Based on pressure and shear stress measurements on the bed, Rajaratnam and Beltaos [8] reported that the most extreme hydrodynamic action is located in the impinging region with maximum pressure at the stagnation point and maximum shear stress just away from it.
On erodible beds, two main crater shapes have been identified [1,2,[7][8][9][10]15]. The first one consists of wide and shallow craters for high impingement heights L (distance from the jet nozzle to the initial flatbed). Here, the flow streamlines weakly deflect in the impinging region, and a smooth crater shape remains stable after the jet flow stop, see Figure 1a. For the second shape, narrower and deeper craters are obtained for small L and large flow rates. Here, the flow streamlines strongly deflect, and the crater shape changes after the jet stop with granular avalanches lowering the steep slopes, see Figure 1b. Here, H and D are the crater depth and width, respectively, and L is the distance above the granular bed.
Aderibigbe and Rajaratnam [9] state that these two crater regimes depend on the value of an adimensional erosion parameter. This parameter is defined for quasi-3D axisymmetric jets as where U J is the mean velocity at the nozzle outlet, b is the nozzle diameter, L is the impingement height, d is the sediment grain size, and s = ρ p /ρ f is the particle/fluid density ratio. This parameter is sometimes called a densimetric Froude number [17] and roughly corresponds to the Rouse number [3]. For the 3D axisymmetric experimental data [9], the erosion parameter E is ranged up to 5. This parameter was used by Sutherland and Dalziel [14] to investigate the crater formation dynamics in granular beds by a 3D water jet. In the 2D turbulent plane jet case, the erosion threshold is well governed by E, as shown in [17]. Recently, Badr et al. [15] proposed that the modified Shields can control the erosion threshold for both turbulent and laminar regimes. However, we should consider the correct decay law of self-similar regimes and the precise position of the corresponding virtual origin. Moreover, in the 2D plane jet test case of Badr et al. [15], the critical value of the local Shields number for erosion is a constant Sh c ≈ 1 independently of the jet flow regime and the grain diameter of the bed.
On the other hand, several numerical studies have been proposed to study the crater formation by an incident jet. We can clearly identify two types of numerical modeling frameworks. In the first type, numerical models apply computational fluid dynamics (CFD) software for computing the fluid motion only above the bed [1,[25][26][27]. The second modeling framework is based on two-phase approaches, in which the motion equations solve both solid and fluid phases [28][29][30][31][32][33][34][35].
The common point of models based on CFD software is that the sediment bed is considered an impermeable wall (no-slip conditions). However, their modeling approaches mainly differ in treating the two-phase interaction between the fluid and the sediment bed. For example, Metzger et al. [1] use CFD++ software to study the steady-state flow above different stationary boundary conditions (or crater geometries) but calculate the porous flow inside the sediment bed by Darcy's law and pressure values. Weidner et al. [25] and Mercier et al. [26] use FLUENT to dynamically adapt the computational grid into the bed's evolving geometry; however, they do not consider the sediment below the bed interface.
For the two-phase formulation, the Eulerian approach usually models the fluid phase, but Eulerian or Lagrangian techniques can be used for the solid phase. A Lagrangian model for the solid phase enables a microscopic description of the particle motion. The flowing (eroded) and jammed (non-eroded) regions in the sediment bed are directly simulated, grain by grain. However, this framework suffers from high computational costs due to the number of solid particles. For examples of this formulation, we refer the work of Kuang et al. [28], and more recently the studies of Benseghier et al. [29,30].
On the other hand, in the Eulerian-Eulerian framework, clear fluid and solid particles are treated as two interpenetrating fluids. The fluid-particle and particle-particle interactions are described at the scale of the continuum media [32][33][34][35]. However, if only a liquid-like behavior of solid particles is modeled to handle the motion of the granular phase for the jet problem, the sediment motion becomes unrealistic, particularly in the jammed region [32,33].
Therefore, solid-like and liquid-like behaviors of sediment are essential and need to be considered for more accurate models. In recent years, few efforts have been performed to include this property in the Eulerian approach; for example, Yuan et al. [33] applied an incipient motion theory in their numerical simulations, Uh Zapata et al. [34] proposed a modified inter-phase momentum transfer formulation, and Wang et al. [35] used a model for the solid-fluid mixture combined with a viscoplastic description of soil. Consequently, this paper proposes a unified formulation for continuum mechanics to simulate vertical-jetinduced erosion with a unique local parameter governing the solid/liquid transition for the sediment. This model is based on a technique initially developed by Greenshields and Weller [36] to handle fluid-structure interactions of flexible and impermeable tubes under fluid flows. This paper's originality relies on the experimental results and new numerical model for submerged vertical plane jet-induced craters on horizontal non-cohesive sediment. The remaining document organization is as follows. Section 2 presents the experimental setup to generate craters, measuring their sizes for different impingement distances and jet velocities. The presentation of the two-phase model, using the two-phase Euler-Euler approach, is proposed in Section 3. We give a complete description of the model improvements, which produce both solid-like and liquid-like behaviors of the sediment. Next, Section 4 shows that the correct local Shields number choice for this governing parameter leads to very close agreement with the experimental results. Finally, we present the discussion, conclusions, and future work.

Laboratory Experiments
In this section, the experimental results for the crater depth and width are analyzed following Badr et al. [15] in terms of the self-similar jet model. We pay particular attention to the corresponding virtual origin. Simple scaling for the crater sizes is found in terms of a local erosion parameter.

Experimental Setup
The experimental set-up consists of a vertical plane water jet impinging a horizontal granular bed immersed in water, as sketched in Figure 2. The grains are monodisperse glass beads with diameter d = 350 µm with the density ρ p = 2.5 × 10 3 kg/m 3 . The granular bed is 10 cm in height (about 300-grain diameters) with a solid volume fraction of about 0.6. It lies on the bottom of a 3D rectangular container of dimensions 50 × 20 × 3.2 cm 3 filled with water (density ρ f = 10 3 kg/m 3 and kinematic viscosity ν f = 10 −6 m 2 /s); see Figure 2a.   A vertical plane water jet is generated at the exit of an injector with an inner rectangular cross-section of dimensions b × w J = 0.4 × 2.6 cm 2 . The outside cross-section of the injector is 1 × 3.2 cm 2 to fit within the container gap. The injector is at a distance L above the granular bed, see Figure 2b.
The mean velocity of the jet is given by U J = Q/(bw J ). A gear pump controls the steady water flow Q with very low negligible fluctuations in a closed circuit. Here, water is sucked through the bed from the porous bottom of the container and re-injected by the pump through the injector.
In this paper, the jet Reynolds number is defined as Although Re is not much relevant to describe the erosion phenomena, the jet Reynolds number is essential to characterize the jet regime among laminar, oscillating laminar, oscillating turbulent, and fully turbulent, as detailed by Badr et al. [15] and Tritton [19]. The injector is long enough (20 cm) to ensure a laminar regime inside the injector. It also has a parabolic Poiseuille velocity profile across the thickness b at the nozzle in the range of jet Reynolds numbers Re < 2 × 10 3 .

Measure of the Crater Dimensions
The erosion process is observed until the resulting crater reaches a quasi-stationary state in a few seconds. The size of the jet-induced scour hole is characterized by the two measured parameters: crater depth and width, as shown in Figure 1. The depth, H, is defined as the distance from the bottom of the crater to the initial bed surface. The width, D, is also determined at the level of the initial bed surface and not at the level of the crater crests.
The measurements are extracted from the images taken by a camera placed in front of the container with backlighting. When the crater stops evolving, we measure the results from the difference between the quasi-stationary image and the initial flatbed image. This procedure allows an easy and accurate determination of H and D. It is important to remark that the measured crater has dynamical sizes, which result from the dynamical crater equilibrium sustained by the steady jet flow. In that case, the dynamic angle may be larger than the avalanche angle.

Experimental Results
The crater size measurements are shown in Table 1 and Figure 3. We use six different jet velocities in the range of 0.16 < U J < 0.47 m/s (630 < Re < 1900) and different jet-bed distances (2 ≤ L/b ≤ 40). Figure 3 shows the crater depth H and width D normalized by the jet thickness b.   Figure 3a shows that for a given mean jet velocity U J , the crater is deeper at lower L values. As one would expect, H decreases toward zero as L increases. Above a critical distance, L c , there is no erosion and thus no crater formation is seen anymore. Note that the critical distance L c also increases with U J . For a given value of L (see for instance L/b = 10), H increases with the jet velocity U J . At low enough U J , here, U J < 0.1 m/s, no erosion can be seen, whatever the value of the jet-bed distance. On the other hand, Figure 3b shows that the crater width D increases with L and U J from a non-zero value, D 0 , close to the erosion threshold.
Previous behavior can be understood when considering the jet spatial evolution downstream of the injector nozzle. As the jet Reynolds number is larger than 200, the jet regime is turbulent with a local velocity u 0 at the jet axis. This velocity decreases according to the distance L with a law that can be approximated by the following self-similar model for free jets [37]: where K ≈ 1/2 is the decay rate and λ is the distance of the virtual origin of the jet from the nozzle exit. Here λ > 0 for a virtual origin downstream of the jet nozzle. In the proposed experimental conditions, λ is equal to 10b. From jet momentum conservation, the jet width δ increases in the meantime with L by the following law: where K 2 ≈ 1/4 corresponds to the jet opening angle θ = tan −1 (K 2 ) ≈ 15 • . We remark that the jet velocity spatial profile (3) and (4) is only valid downstream the potential core. By contrast, in the injector, the jet velocity u 0 remains constant and equal to u 0 = (3/2)U J for a laminar parabolic Poiseuille profile. In such a self-similar model, the reduced jet-bed distance, L − λ, appears as the pertinent length scale that may govern the crater size. Thus, taking into account the virtual origin of the jet, a relevant dimensionless parameter for the jet erosion strength is the local erosion parameter E or Rouse number. In the present inertial regime of high particle Reynolds number, this parameter is given by which corresponds to the ratio of the local jet velocity to the grain settling velocity. Equation (5) is similar to expression (1) proposed in [9]. Note that if U J goes with √ b in the numerator of Equation (5), then non-dimensionalization will need a √ L factor in the denominator. It is all a direct consequence of the jet being characterized by its momentum flux. Thus, besides considering the effective jet-sediment distance, the model differs in the exponent −1/2 (instead of −1) for the L scaling. This exponent is coming from the self-similar free jet models for a 2D plane (instead of 3D axisymmetric) flow geometry. Figure 4 shows the crater size data normalized by the reduced distance L − λ as a function of the dimensionless erosion parameter E given by Equation (5). We compute E using s − 1 = 1.5, g = 9.8 m/s 2 and d = 350 µm. Note that the different curves collapse into a master curve of simple trend. The dimensionless crater depth H/(L − λ) increases linearly with E above a critical value E c corresponding to erosion threshold, as follows: where m H = 0.7 ± 0.05 and E c = 1 ± 0.1. Similarly, the dimensionless crater size D/(L − λ) increases linearly with E above the critical value E c ; however, it has a non-zero value D 0 at the threshold. Thus, we obtain where The critical value E c ≈ 1 found here for the erosion threshold at vanishing crater depth agrees very well with the critical value one reported by Badr et al. [15] with just a visual criterion of the few first moving grains in a similar setup. It is worth noting that taking into account the distance λ of the virtual origin of the jet from the nozzle exit is crucial to obtain a good rescaling of the data and such simple scaling laws. Indeed, the critical value E c ≈ 1 at the threshold is very satisfying, meaning that E is a pertinent parameter for describing vertical jet-induced erosion.
Finally, considering the evolution of crater sizes H and D with the jet flow action, the shape of the crater evolves with a typical slope H/(D/2) that tends toward the value 2m H /m D ≈ 1 at large erosive strength E (E 1). This value corresponds to a global slope angle of about 45 • , thus much larger than the usual critical avalanche angle 25 • without any jet flow [38]. This shows that at large E, the crater slope is enhanced by an upward-back flow of the jet, as happens in homogeneous parallel shear flow [39].

Mathematical Model
The present two-phase flow model comes from a strict application of the mathematical formulation of Drew and Lahey [40], using a Eulerian-Eulerian description for the fluidparticle system. In this section, the governing equations for the two phases corresponding to clear fluid (water) and disperse medium (sediment particles in a fluid) are briefly recalled with special closure laws.
The present two-phase flow model is initially based on the two-fluid model presented in Barbry et al. [41] and Nguyen et al. [42,43]. This model only assumes a liquid-like behavior for the disperse phase (solid grains). We refer to this formulation as the standard or classical model. However, it cannot account for a possible solid-like behavior of the granular bed. In the present study on scouring erosion, we propose an improved model considering mobile and immobile parts of the granular media by a simple switch from fluid mechanics to soil mechanics equations. The benefit of this new treatment for the sediment particles is illustrated and justified in the numerical results when the method is compared with the initial two-fluid model.

Governing Equations
Before specifying the novel feature introduced to the model in the present paper, let us first recall the initial two-fluid model. For a constant width-integrated 2D x/z two-phase flow model, the governing equations in a Eulerian formulation for the k-phase read as where k = f stands for the clear fluid (water) and k = s for the solid particles, respectively. The volume fraction is represented by values α k with α f + α s = 1. Here, ρ k is the density, u k = (u k , w k ) is the velocity vector, and g stands for gravity. The symbol M k refers to the inter-phase momentum transfer, andσ k is the total stress tensor of phase k. In this paper, the double bar over a symbol means a tensor of order two. The total stress tensor of phase k is given bȳ whereĪ is the identity tensor. Equation (10) is the sum of three contributions: the pressure (p k ) of phase k, the deviatoric viscous (τ k ), and Reynolds stresses (τ Re k ). Following Lundgren [44] for the slow motion of two-fluids, the non-turbulent stresses are expressed by whereD k = 1 2 ∇u k + ∇u T k is the shear rate tensor of phase k. In Equations (11) and (12), the parameter µ f is the dynamic viscosity of the clear fluid (water) and β is an increasing function of the solid volume fraction of grains, for which we use the expression of Graham [45]. We refer to Nguyen et al. [42,43] for more details.
To compute the Reynolds stressesτ Re k for the two-phase model, the standard k f − f model is used for the fluid phase. For the turbulent model of the solid phase, we use the closure equations k s − k f s (turbulent kinetic energy of particulate phase, k s , with covariance of fluid and solid velocity fluctuations, k f p ). This model is considered a twoway coupling because it can evaluate the effects of solid particles on fluid turbulence and vice versa [42,46].
The momentum exchange between phases, M k , is decomposed into jump conditions (p ki andτ ki ) and different forces (M * k ) acting on phase k, as follows Here, the subscript i specifies the solid/fluid interface. As the interfacial tension is zero, there is no discontinuity of the pressure (p si = p f i = p f − ρ f ||u f − u s || 2 /4, where · means the Euclidean norm) and the shear stress (τ si =τ f i = βτ f ) across the fluid/sediment interface (Drew and Lahey [40]).
In Equations (13) and (14), the momentum exchanged between phases (M * f = −M * s ) is the sum of the drag and lift forces exerted on the solid particles, the added mass effect, the Faxén correction and the Basset historical term. In the present formulation, we only retain the drag force as it is dominant, as follows: where u r = u f − u s − u d is the relative velocity of fluid-particles. Here, u d is the drift velocity which represents the correlation between the fluctuating velocity of the fluid and the instantaneous spatial distribution of particles. The drag force is also expressed as a function of the solid concentration (α s ) and the particle relaxation time τ f s given by For turbulent two-phase flows, the drift velocity and particle relaxation time are very important, as discussed in Nguyen et al. [42]. For the present study, the drag coefficient (C D ) is computed from the formula of Haider and Levenspiel [47], where Re p is the particle Reynolds number.

A Unified Momentum Equation for the Solid Phase
The set of Equations (8)-(16) roughly describes the fluid-sediment system as a twofluid system. In particular, Equations (11) and (12) for the viscous stress correspond to a rheological law for complex fluids. Such a model is only valid below the packing fraction when the granular material has a liquid-like suspension behavior. However, we expect a solid-like behavior for the particulate phase close to the packing fraction at the bed surface or deeper in bulk.
The numerical solution given by Equations (8)- (16) are unphysical and erroneous in high solid fraction regions: the sediment's liquid-like behavior is clearly visualized by wavy iso-contour lines in the sediment bed; see Section 4.1. We remark that similar results are obtained using the Eulerian model proposed by Qian et al. [32] and Yuan et al. [33] for this problem. To overcome this issue, this paper adapts a unified formulation, mainly inspired by Greenshields and Weller [36], to describe both the liquid-like and solid-like behavior of the sediment phase.

Solid-Like Model
For the solid-like model, we assume a Hookean material that undergoes small deformations. In this case, a linear elastic, isotropic model is used. Thus, the strain tensor (ε s ) corresponds to the time integration of the strain rate tensor (D s ) as follows: where the deviatoric component ofσ is defined as dev(σ) ≡σ − 1 3 tr(σ)Ī, and dev(σ 0 ) is the deviatoric solid stress at the initial time. Equation (17) can be approximated by a quadrature formula for the integral, as follows: where ω i represents the quadrature weighting functions. Here, (·) n indicates the current time (t n ). Note that the sum is taken over previous stages t i (i = 1, . . . , n − 1). Next, we can substitute (18) into (17). Now, using the relationship between ε s and D s , we obtainσ where the shear modulus is given in terms of the elastic Young modulus (E) and the Poisson's ratio (ν) by G = E/2(1 + ν). In Equation (19), Σ is given by representing the accumulation of elastic stress from the beginning of the loading. Thus, the solid-like behavior of the granular phase can be modeled as Finally, we obtain a unified momentum equation for the sediment incorporating (21) into the two-phase flow model (8)- (15). Thus, to unify liquid-like model (12) and solid-like model (21), we introduce a smooth transition function (F), as follows: where 0 ≤ F ≤ 1. We remark that if F = 1 then Equation (22) reduces to the solid-like model (21). If F = 0 then a liquid-like motion is obtained by Equation (12). In other terms, the improvement proposed in this paper has a non-intrusive character since the original two-fluid model is recovered. For the fluid phase, the equations remain unchanged.

Solid-Liquid Transition of the Granular Phase
In the original development of Greenshields and Weller [36], F is a step function that takes either the value 0 or 1 according to the fluid or solid regions, respectively, since there is no interpenetration between them. In our modeling, F varies continuously from 1 to 0 using a smooth solid-liquid transition for the sediment.
From a physical viewpoint, function F encodes the erosion or fluidization criterion because it determines at any point of the computational domain if the disperse phase is flowing (eroded/fluidized bed for F = 0) or not (static bed for F = 1). Therefore, this function's formulation is as complex as the physics of the erosion threshold. From a numerical point of view, function F should be monotone and sufficiently smooth with a continuous first derivative to avoid any numerical instability. This paper considers for F an S-shape function with the following analytical expression: where erf is the Gaussian error function. We remark that Equation (23) has a unique parameter ξ governing the liquid-solid transition. In the numerical simulations, we select corresponding to the deviation of the local Shields number Sh from its critical value for erosion Sh c ≈ 1. This choice is consistent with the erosion criteria of Badr et al. [15] for the plane jet impinging into granular beds. This selection also agrees with the dimensionless crater depth and width formulas given in the experimental results, as presented in the following section.

Numerical Technique
To solve the system of Equations (8)- (17), we use the techniques initially developed by Guillou and Barbry [41,48] and later improved by Uh Zapata et al. [34,49]. This in-house program is currently written in Fortran 90, displaying the results in Paraview software.
The σ-transformation is applied to the vertical coordinates to fit the computing mesh to the free surface evolution. A finite volume formulation approximates spatial derivatives. A projection technique [50] is applied to calculate the pressure and velocity for each phase. We use a staggered grid to avoid the spurious oscillations induced by the projection technique [51]. Advection terms in Equation (9) are handled by a total variation diminishing (TVD) scheme. The time scheme is implicit in both vertical and horizontal directions.
A successive order relaxation (SOR) iterative solver is used for the resulting linear equation system. Furthermore, the code was fully parallelized through a message passing interface (MPI) and graphics processing unit (GPU), which significantly optimize the performance of the two-phase flow program.

Numerical Results
In this section, the numerical results are obtained using a 2D rectangular domain with 0.25 m of height and 0.2 m of width. We impose a Poiseuille profile of velocity at the nozzle outlet of thickness b, which is fixed at a distance L from the bed.
In the simulations, we consider a rectangular mesh grid. The highest resolution is given by a uniform mesh of 401 points along the horizontal and 501 points in the vertical direction. Thus, the spatial steps are given by ∆x = ∆z = 0.5 mm. This choice corresponds to a grid resolution close to one-grain diameter. It would not be pertinent to consider finer meshes as we can not have grid resolutions more refined than the sediment diameter for the Euler-Euler two-phase flow formulation. For this fine space resolution, small time steps are reduced to ∆t = 2 × 10 −6 s to have stable simulations. In agreement with the experiments, each run lasts a few seconds of simulated time.

Original and Unified Model
To illustrate the benefit of the proposed treatment for the solid-liquid transition, Figure 5 presents the numerical results using the original two-phase flow model (transition function F = 0). For these simulations, we use a coarse grid (∆x = 1 mm, ∆z = 2 mm), and parameters b = 4 mm, U J = 0.4 m/s, and L = 2 cm). The initial flatbed interface is located at z = −0.15 m. Although we have a crater at the early stages, the bed shows unrealistic wavy modulation of the bed.
On the other hand, Figure 6 shows the numerical results using the proposed unified two-phase Euler-Euler model. For these simulations, we use the same parameters as the standard model. Note that the qualitative behavior is now close to the experimental observations: the bed does not present a wavy behavior but a localized crater. Furthermore, zero particle velocity is predicted within the sediment bed. In other words, our numerical treatment for the solid-liquid transition in the model mainly corrects the identified numerical inaccuracies associated with the original two-fluid formulation.
As previously pointed out, the profile function plays an essential role in modeling correctly the solid phase. Figure 7 shows the profiles for F given by Equation (23) and the numerical results of previous simulations. Note how this function correctly varies continuously from 1 to 0 dividing the computational domain into liquid-like and solid-like regions. We remark the correlation between the scour dimensions and the values of F. Although previous results shows the correct implementation of the model, the mesh resolution does not allow to fully observe the complex structure of the water jet and crater formation. Figures 8 and 9 present the velocity field and solid volume fraction results at t = 1.5 seconds using b = 8 mm, U J = 0.323 m/s, and L = 4.6 cm. In these simulations, we use a fine (∆x = ∆z = 0.5 mm), a medium (∆x = ∆z = 1 mm), and a coarse (∆x = ∆z = 2 mm) grid resolutions. As illustrated in Figure 8, the plane jet flow shows a complex structure with some fluctuation and vortices along the jet. However, the jet flow structure is out of the scope of the present paper, so we only focus on the cratering. Note that the region of solid-like behavior of the sediment is characterized by low fluid velocity, in contrast to the high velocity just above this region.
The abrupt variation in the velocity magnitude indicates that a continuous interface line can model the solid-liquid transition zone. This interface cane be verified in Figure 9 showing 100 contour lines of the solid volume fraction distributed uniformly in the interval [0, 0.55]. Note that the transition zone between the sediment bed and the liquid is less diffusive as the mesh resolution increases, as expected. This interface is used to measure the crater size (H, D).

Crater Size Predictions
As discussed in Section 3.3, we use the local parameter ξ = Sh − Sh c for governing the sediment solid-liquid transition, where Sh is the local Shields number Sh and Sh c stands for its critical value. Notice that this formulation considers the total fluid velocity instead of only its tangential component. Indeed, we systematically obtain a zero value of the local Shields parameter at the stagnation point in front of the jet by considering only the horizontal component of the fluid velocity. As a result, the scour hole geometry would present a non-eroded zone at the jet axis with a whiskers-like shape.
Numerical results have shown to have an excellent agreement with experimental data using Sh c = 0.9. This value also agrees with the critical Shields number near one reported by Badr et al. [15]. Figure 10 presents the fluid velocity field for b = 4 mm with different values of the impingement distance (L) and mean jet velocity (U J ) using the same Sh c = 0.9. Results demonstrate that the crater is wider and shallower as L increases, as expected. We also notice that the velocity field for the fluid phase presents more oscillations for higher jet velocities. Experiments also show these jet oscillations. Moreover, the proposed two-phase model generates both strongly and weakly deflected cases with recirculation zones inside or outside the crater, respectively (see Figure 10c,d). To verify the assumptions made for the velocity based on the Shield parameter and its critical value, we test a series of three numerical simulations for different jet widths b equal to 4 mm, 6 mm, and 8 mm ( Table 2). The critical value Sh c = 0.9 is kept constant for all numerical simulations. Numerical results confirm that D is larger and H is smaller as L increases. We note that the crater dimensions also increases as the mean jet velocity U J is larger, as expected. Finally, Figure 11 compares the numerical and experimental results of the crater dimensions D and H. The experimental data are the same as those presented in Figure 4. The straight lines correspond to the linear models from experimental data given by Equations (6) and (7). The numerical results are made according to the information given in Table 2. The comparison shows a very good agreement for all values. These findings show the ability of the proposed two-phase model to reproduce complex flows with sediment. Figure 11. Comparison between experiments and numerical results for (a) the crater depth and (b) width normalized by the reduced jet-bed distance (L − λ). Solid lines correspond to linear laws given by Equations (6) and (7).

Discussion and Conclusions
In this paper, we study experimentally and numerically the crater sizes generated on horizontal non-cohesive sediment by a submerged vertical 2D plane jet. We showed that the crater depth and width, measured under various impinging distances and jet velocities conditions, are governed by the effective impinging distance L − λ corrected from the virtual jet origin, and by a dimensionless erosion parameter E. This erosion parameter, inspired by the one proposed by Aderibigbe and Rajaratnam [9], also accounts for the virtual origin of the jet, as proposed by Badr et al. [15] and Sutherland and Dalziel [14]. The dimensionless crater sizes evolve approximately linearly with respect to E above the critical value E c ≈ 1 corresponding to the erosion threshold. In contrast to the crater depth, the crater width is non-zero at the threshold, as already found in [10,14].
A two-phase numerical model was used successfully to simulate the crater formation. We suppress numerical inaccuracies from a standard two-fluid model by considering a new formulation for the sediment. It consists of a unified equation for the solid phase, which employs a smooth function governing the transition between the liquid-like and solid-like regimes. This function is based on the deviation of a local Shields number from a critical value Sh c ≈ 1. Numerical results agree accurately with the experimental data for all the proposed tests. Further studies will focus on testing the proposed approach with different configurations. For instance, we are interested in scour holes generated by inclined or even horizontal jet flows.

Acknowledgments:
The authors extend special thanks to Compute Canada (contract No. 3148) for providing access to the computing facility.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: 2D Two