Damien Pham-Van-Bang

Miguel Uh Zapata

Georges Gauthier

Philippe Gondret

Wei Zhang

Dan Kim

Nguyen

Damien Pham Van Bang

Kim Dan Nguyen

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

Keywords: jet erosion test, sediment transport, numerical modeling, two-phase flow model, experiments, water injection dredging

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Vertical jet-induced scour erosion of soil has been studied for many industrial applications such as aerospace or hydraulic engineering [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Metzger | Jet-induced cratering of a granular surface with application to lunar spaceports[END_REF][START_REF] Rouse | Criteria for similarity in the transportation of sediment[END_REF][START_REF] Hanson | Apparatus, test procedures, and analytical methods to measure soil erodibility in situ[END_REF][START_REF] Hanson | Lessons learned using laboratory JET method to measure soil erobibility of compacted soils[END_REF][START_REF] Perng | Underwater sand bed erosion and internal jump formation by a travelling plane jets[END_REF]. As reported by Metzger et al. [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Metzger | Jet-induced cratering of a granular surface with application to lunar spaceports[END_REF], 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 [START_REF] Rouse | Criteria for similarity in the transportation of sediment[END_REF] for testing criteria on erosion, the study of Hanson and Cook [START_REF] Hanson | Apparatus, test procedures, and analytical methods to measure soil erodibility in situ[END_REF] assessing in situ the erodibility of soil material, and the research developed by Hanson and Hunt [START_REF] Hanson | Lessons learned using laboratory JET method to measure soil erobibility of compacted soils[END_REF] in the so-called Jet Erosion Test (JET). We also have the work of Perng and Capart [START_REF] Perng | Underwater sand bed erosion and internal jump formation by a travelling plane jets[END_REF] for Water Injection Dredging (WID) and trenching in harbors by jet translation and possible jet inclination.

Experimental investigations of the jet normal impingement on plane surfaces have been conducted in various conditions. Although the jet nozzle section is circular in most cases [START_REF] Hanson | Apparatus, test procedures, and analytical methods to measure soil erodibility in situ[END_REF][START_REF] Kobus | Flow field and scouring effects of steady and pulsating jets impinging on a movable bed[END_REF][START_REF] Rajaratnam | Erosion by impinging circular turbulent jets[END_REF][START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF][START_REF] Rajaratnam | Erosion of Sand by Circular Impinging Water Jets with Small Tailwater[END_REF][START_REF] Rajaratnam | Impingement of circular turbulent jets on rough boundaries Impact des jets turbulents sur des parois rugueuses[END_REF][START_REF] Mazurek | Assessment of the Erodibility of a Cohesive Soil using a Submerged Circular Turbulent Impinging Jet[END_REF][START_REF] Sutherland | Bedload transport by a vertical jet impinging upon sediments[END_REF], it can also be rectangular [START_REF] Beltaos | Plane turbulent impinging jets[END_REF][START_REF] Rajaratnam | Erosion by plane turbulent jets[END_REF][START_REF] Phares | The wall shear stress produced by a normal impingement of a jet on a flat surface[END_REF][START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF]. The circular and rectangular cases are called quasi-three-dimensional (3D) axisymmetric and two-dimensional (2D) plane jets, respectively. The jet can be either gas [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Metzger | Jet-induced cratering of a granular surface with application to lunar spaceports[END_REF] or liquid [START_REF] Hanson | Apparatus, test procedures, and analytical methods to measure soil erodibility in situ[END_REF][START_REF] Perng | Underwater sand bed erosion and internal jump formation by a travelling plane jets[END_REF] within unsubmerged [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Metzger | Jet-induced cratering of a granular surface with application to lunar spaceports[END_REF] or submerged [START_REF] Rajaratnam | Erosion by impinging circular turbulent jets[END_REF][START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF][START_REF] Rajaratnam | Erosion of Sand by Circular Impinging Water Jets with Small Tailwater[END_REF][START_REF] Rajaratnam | Impingement of circular turbulent jets on rough boundaries Impact des jets turbulents sur des parois rugueuses[END_REF][START_REF] Mazurek | Assessment of the Erodibility of a Cohesive Soil using a Submerged Circular Turbulent Impinging Jet[END_REF][START_REF] Sutherland | Bedload transport by a vertical jet impinging upon sediments[END_REF][START_REF] Beltaos | Plane turbulent impinging jets[END_REF][START_REF] Rajaratnam | Erosion by plane turbulent jets[END_REF][START_REF] Phares | The wall shear stress produced by a normal impingement of a jet on a flat surface[END_REF][START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF] conditions. Finally, the soil can be either cohesive [START_REF] Rajaratnam | Erosion of Sand by Circular Impinging Water Jets with Small Tailwater[END_REF] or non-cohesive [START_REF] Mazurek | Assessment of the Erodibility of a Cohesive Soil using a Submerged Circular Turbulent Impinging Jet[END_REF][START_REF] Beltaos | Plane turbulent impinging jets[END_REF].

Based on experiments of 3D axisymmetric jets impinging on a smooth non-erodible wall, Rajaratnam and Beltaos [START_REF] Rajaratnam | Erosion by impinging circular turbulent jets[END_REF] and Phares et al. [START_REF] Phares | The wall shear stress produced by a normal impingement of a jet on a flat surface[END_REF] divide the flow into three zones:

•

The free jet region. There is no influence on the boundary. It can be described by self-similar solutions [START_REF] Tritton | Physical fluid dynamics[END_REF][START_REF] Pope | Turbulent flows[END_REF].

• The impinging region. It is characterized by a stagnation point on the jet axis at the wall with a deflexion of flow streamlines [START_REF] Donalson | A study of free jet impingement. Part 1. Mean properties of free and impinging jets[END_REF].

•

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 [START_REF] Glauert | The wall jet[END_REF][START_REF] Launder | The turbulent wall: Measurements and modeling[END_REF][START_REF] Magyari | The algebraically decaying wall jet[END_REF].

Based on pressure and shear stress measurements on the bed, Rajaratnam and Beltaos [START_REF] Rajaratnam | Erosion by impinging circular turbulent jets[END_REF] 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 [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Metzger | Jet-induced cratering of a granular surface with application to lunar spaceports[END_REF][START_REF] Kobus | Flow field and scouring effects of steady and pulsating jets impinging on a movable bed[END_REF][START_REF] Rajaratnam | Erosion by impinging circular turbulent jets[END_REF][START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF]. 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. 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. Aderibigbe and Rajaratnam [START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF] 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

E = U J b L 1 (s -1)gd 1/2 , ( 1
)
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 [START_REF] Rajaratnam | Erosion by plane turbulent jets[END_REF] and roughly corresponds to the Rouse number [START_REF] Rouse | Criteria for similarity in the transportation of sediment[END_REF].

For the 3D axisymmetric experimental data [START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF], the erosion parameter E ranged up to 5. This parameter has been used by Sutherland and Dalziel [START_REF] Sutherland | Bedload transport by a vertical jet impinging upon sediments[END_REF] 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 [START_REF] Rajaratnam | Erosion by plane turbulent jets[END_REF]. Recently, Badr et al. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF] 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. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF], 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 [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF][START_REF] Weidner | Numerical simulation of jet test and associated soil erosion[END_REF][START_REF] Mercier | On the numerical modelling of the Jet Erosion Test[END_REF][START_REF] Wang | Study on deepwater conductor jet excavation mechanism in cohesive soil[END_REF]. The second modeling framework is based on two-phase approaches, in which the motion equations solve both solid and fluid phases [START_REF] Kuang | Discrete particle simulation of jet-induced cratering of a granular bed[END_REF][START_REF] Benseghier | Relevance of free jet model for soil erosion by impinging jets[END_REF][START_REF] Benseghier | A parallel GPU-based computational framework for the micromechanical analysis of geotechnical and erosion problems[END_REF][START_REF] Boyaval | A semi-Lagrangian splitting method for the numerical simulation of sediment transport with free surface flows[END_REF][START_REF] Qian | Numerical simulation of sediment erosion by submerged jets using an Eulerian model[END_REF][START_REF] Yuan | Numerical study of sand scour with a modified Eulerian model based on incipient motion theory[END_REF][START_REF] Zapata | Parallel simulations for a 2D x/z two-phase flow fluid-solid particle model[END_REF][START_REF] Wang | Modeling the hydraulic excavation of cohesive soil by a moving vertical jet[END_REF].

The common point of models based on CFD software is that the sediment bed is considered as impermeable wall (no-slip conditions). But their modeling approaches mainly differ in treating the two-phase interaction between the fluid and the sediment bed. For example, Metzger et al. [START_REF] Metzger | Craters formed in granular beds by impinging jets of gas[END_REF] 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. [START_REF] Weidner | Numerical simulation of jet test and associated soil erosion[END_REF] and Mercier et al. [START_REF] Mercier | On the numerical modelling of the Jet Erosion Test[END_REF] 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. [START_REF] Kuang | Discrete particle simulation of jet-induced cratering of a granular bed[END_REF], and more recently the studies of Benseghier et al. [START_REF] Benseghier | Relevance of free jet model for soil erosion by impinging jets[END_REF][START_REF] Benseghier | A parallel GPU-based computational framework for the micromechanical analysis of geotechnical and erosion problems[END_REF].

On the other hand, in the Eulerian-Eulerian framework, clear fluid and solid particles are treated as two interpenetrating fluids. The fluid-particles and particle-particle interactions are described at the scale of the continuum media [START_REF] Qian | Numerical simulation of sediment erosion by submerged jets using an Eulerian model[END_REF][START_REF] Yuan | Numerical study of sand scour with a modified Eulerian model based on incipient motion theory[END_REF][START_REF] Zapata | Parallel simulations for a 2D x/z two-phase flow fluid-solid particle model[END_REF][START_REF] Wang | Modeling the hydraulic excavation of cohesive soil by a moving vertical jet[END_REF]. 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 [START_REF] Qian | Numerical simulation of sediment erosion by submerged jets using an Eulerian model[END_REF][START_REF] Yuan | Numerical study of sand scour with a modified Eulerian model based on incipient motion theory[END_REF]. 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 has been performed to include this property in the Eulerian approach; for example, Yuan et al. [START_REF] Yuan | Numerical study of sand scour with a modified Eulerian model based on incipient motion theory[END_REF] apply an incipient motion theory in their numerical simulations, Uh Zapata et al. [START_REF] Zapata | Parallel simulations for a 2D x/z two-phase flow fluid-solid particle model[END_REF] proposed a modified inter-phase momentum transfer formulation, and Wang et al. [START_REF] Wang | Modeling the hydraulic excavation of cohesive soil by a moving vertical jet[END_REF] use 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-jet induced 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 [START_REF] Greenshields | A unified formulation for continuum mechanics applied to fluid-structure interaction in flexible tubes[END_REF] 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. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF] 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 2(a). The container walls are in PolyMethyl MethAcrylate (PMMA), allowing visualization.

A vertical plane water jet is generated at the exit of an injector of 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 2(b).

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

Re = U J b ν f . ( 2
)
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. 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, as shown in Figure 1.

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 results 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 andFigure 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.

Table 1. Experimental results using different experimental parameters of U J , L and jet thickness b = 0.004 m. Figure 3(a) shows that for a given mean jet velocity U J , the crater depth H is higher at lower L values. Furthermore, H decreases towards 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 3(a) shows that the crater width D increases with L and U J for 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 [START_REF] Chua | Measurements of a confined jet[END_REF]:

Mean velocity

Reynolds number

Impingement distance

Crater width

Crater depth

Mean velocity

Reynolds number

Impingement distance

Crater width

Crater

depth U J (m/s) Re L (m) D (m) H (m) U J (m/s) Re L (m) D (m) H (
u 0 U J = 1 K √ L/b -λ/b , (3)
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:

δ b = K 2 (L/b -λ/b), ( 4
)
where K 2 ≈ 1/4 corresponds to the jet opening angle θ = tan -1 (K 2 ) ≈ 15 o . We remark that the jet velocity spatial profile (3)-( 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

E = U J b L -λ 1/2 1 (s -1)gd 1/2 , ( 5
)
which corresponds to the ratio of the local jet velocity to the grain settling velocity. Equation ( 5) is similar to expression (1) proposed in [START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF]. However, 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 [START_REF] Hanson | Lessons learned using laboratory JET method to measure soil erobibility of compacted soils[END_REF]. 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

H/(L -λ) = m H (E -E c ), (6)
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 threshold. Thus we get

D/(L -λ) = m D (E -E c ) + D * 0 , (7)
where m D = 1.4 ± 0.1 and D * 0 = D 0 /(Lλ) = 0.55 ± 0.05. 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. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF] 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 o , thus much larger than the usual critical avalanche angle 25 o without any jet flow [START_REF] Midi | On dense granular flows[END_REF]. 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 [START_REF] Loiseleux | Onset of erosion and avalanche for an inclined granular bed sheared by a continuous laminar flow[END_REF].

Mathematical Model

The present two-phase flow model comes from a strict application of the mathematical formulation of Drew and Lahey [START_REF] Drew | Analytical Modelling of Multiphase Flow[END_REF] using an 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. [START_REF] Barbry | Une approche diphasique pour le calcul du transport sédimentaire en milieux estuariens[END_REF] and Nguyen et al. [START_REF] Nguyen | A two-phase numerical model for suspended-sediment transport in estuaries[END_REF][START_REF] Nguyen | Simulation of dredged sediment releases into homogeneous water using a two-phase model[END_REF]. 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 an Eulerian formulation for the k-phase read:

∂(α k ρ k ) ∂t + ∇ • (α k ρ k u k ) = 0, (8)
∂(α k ρ k u k ) ∂t + ∇ • (α k ρ k u k u k ) = ∇ • (α k σk ) + α k ρ k g + M k , (9)
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 by σk

= -p k Ī + τk + τRe k . ( 10
)
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 [START_REF] Lundgren | Slow flow through stationary random beds and suspensions of spheres[END_REF] for the slow motion of two-fluids, the non-turbulent stresses are expressed by:

α f σf = -α f p f Ī + 2µ f α s Ds + α f D f , (11)
α s σs = -α s p s Ī + 2µ f α s β α s Ds + α f D f , (12)
where Dk = 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 [START_REF] Graham | On the viscosity of suspensions of solid spheres[END_REF]. We refer to Nguyen et al. [START_REF] Nguyen | A two-phase numerical model for suspended-sediment transport in estuaries[END_REF][START_REF] Nguyen | Simulation of dredged sediment releases into homogeneous water using a two-phase model[END_REF] 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 sk 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 two-way coupling because it can evaluate the effects of solid particles on fluid turbulence and vice-versa [START_REF] Nguyen | A two-phase numerical model for suspended-sediment transport in estuaries[END_REF][START_REF] Chauchat | On turbulence closures for two-phase sediment-laden flow models[END_REF].

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

M f = p f i Ī -τf i ∇α f + M * f , (13)
M s = p si Ī -τsi ∇α s + M * s . (14)
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 Euclidean norm) and the shear stress ( τsi = τf i = β τf ) across the fluid/sediment interface (Drew and Lahey [START_REF] Drew | Analytical Modelling of Multiphase Flow[END_REF]).

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

M * s = F D u r , F D = α s ρ s /τ f s , (15)
where u r = u fu su 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

τ f s = 4dρ s / 3ρ f C D ||u r || , C D = 24 Re p f (Re p ). ( 16
)
For turbulent two-phase flows, the drift velocity and particle relaxation time are very important, as discussed in Nguyen et al. [START_REF] Nguyen | A two-phase numerical model for suspended-sediment transport in estuaries[END_REF]. For the present study, the drag coefficient (C D ) is computed from the formula of Haider and Levenspiel [START_REF] Haider | Drag coefficient and terminal velocity of spherical and non-spherical particles[END_REF], 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. [START_REF] Qian | Numerical simulation of sediment erosion by submerged jets using an Eulerian model[END_REF] and Yuan et al. [START_REF] Yuan | Numerical study of sand scour with a modified Eulerian model based on incipient motion theory[END_REF] for this problem. To overcome this issue, this paper adapts a unified formulation, mainly inspired by Greenshields and Weller [START_REF] Greenshields | A unified formulation for continuum mechanics applied to fluid-structure interaction in flexible tubes[END_REF], 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 ( Ds ) as follows:

σ s = -p s Ī + dev( σs ) = -p s Ī + t n t 0 dev( σs )dt + dev( σ0 ), ( 17
)
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

t n t 0 dev(σ s )dt = ω n dev( σs ) n + n-1 ∑ i=1 ω i dev( σs ) i , ( 18
)
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 [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF]. Now, using the relationship between ε s and D s , we get σs = -p s Ī + 2Gω n dev( Ds ) + Σ, [START_REF] Pope | Turbulent flows[END_REF] 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

Σ = n-1 ∑ i=1 ω i dev( σs ) i + dev( σ0 ), (20)
representing the accumulation of elastic stress from the beginning of the loading. Thus, the solid-like behavior of the granular phase can be model as

α s σs = -α s p s Ī + α s 2Gω n dev( Ds ) + α s Σ. (21)
Finally, we obtain a unified momentum equation for the sediment incorporating ( 22) into the two-phase flow model ( 8)- [START_REF] Rajaratnam | Erosion by plane turbulent jets[END_REF]. Thus, to unify liquid-like model ( 12) and solid-like model ( 21), we introduce a smooth transition function (F), as follows

α s σs = (1 -F) -α s p s Ī + α s βµ α s Ds + α f D f (22) +F -α s p s Ī + α s 2Gω n dev( Ds ) + α s Σ ,
where 0 ≤ F ≤ 1. We remark that if F = 1 then equation ( 22) reduces to the solid-like model [START_REF] Glauert | The wall jet[END_REF]. If F = 0 then a liquid-like motion is obtained by equation [START_REF] Mazurek | Assessment of the Erodibility of a Cohesive Soil using a Submerged Circular Turbulent Impinging Jet[END_REF]. In other terms, the improvement proposed in this paper has non-intrusive character since the original two-fluid model is recovered. For the fluid phase, equations remain unchanged.

Solid-Liquid Transition of the Granular Phase

In the original development of Greenshields and Weller [START_REF] Greenshields | A unified formulation for continuum mechanics applied to fluid-structure interaction in flexible tubes[END_REF], 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 a S-shape function with the following analytical expression

F = 1 2 (1 -erf(ξ)), ( 23
)
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

ξ = Sh -Sh c , Sh = u 2 f + w 2 f (s -1)gd . ( 24
)
corresponding to the deviation of the local Shields number Sh from its critical value for erosion Sh c . This choice is consistent with the erosion criteria of Badr et al. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF] 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 et al. [START_REF] Guillou | Calcul numérique des ondes de surface par une méthode de projection et un maillage eulérien adaptatif[END_REF] and later improved by Uh Zapata et al. [START_REF] Zapata | Two-phase experimental and numerical studies on scouring at the toe of vertical seawall[END_REF]. 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 [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] 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 [START_REF] Guillou | An improved technique for solving two-dimensional shallow water problems[END_REF]. 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 has been fully parallelized through 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 fine grid resolution close to one-grain diameter. Thus, it would not be pertinent to consider finer meshes. For this fine space resolution, small time steps are reduced to ∆t = 2 × 10 -6 sec 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 and9 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 Figures 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. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF]. 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 Figures 10(c 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 is the same as the one presented in Figure 4. The straight lines correspond to the linear models from experimental data given by equations ( 6) and [START_REF] Kobus | Flow field and scouring effects of steady and pulsating jets impinging on a movable bed[END_REF]. The numerical results are taking according to the information given in Table 2. The comparison shows a very good agreement for all values. These findings show the capacity of the proposed two-phase model to reproduce complex flows with sediment.

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 have shown that the crater depth and width, measured under various impinging distances and jet velocities conditions, are governed by the effective impinging distance Lλ. This distance is corrected from the virtual jet origin and by a dimensionless erosion parameter E. The proposed erosion parameter, inspired by the one proposed by Aderibigbe and Rajaratnam [START_REF] Aderibigbe | Erosion of loose beds by submerged circular impinging vertical turbulent jets[END_REF], also accounts for the virtual origin of the jet, as proposed by Sutherland and Dalziel [START_REF] Sutherland | Bedload transport by a vertical jet impinging upon sediments[END_REF] and Badr et al. [START_REF] Badr | Erosion threshold of a liquid immersed granular bed by an impinging plane liquid jet[END_REF]. The dimensionless crater sizes evolve approximately linearly with respect to E above the critical value for the erosion threshold (E c ≈ 1). By contrast to the crater depth, the crater width is non-zero at the threshold, as already found in [START_REF] Sutherland | Bedload transport by a vertical jet impinging upon sediments[END_REF].

A two-phase numerical model has been 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 liquidlike 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 to the experimental data for all the proposed tests. Further studies are focused on testing the proposed approach with different configurations. For instance, we are interested in scour holes generated by inclined or even horizontal jet flows.

Figure 1 .

 1 Figure 1. Geometrical parameters of the eroded sediment bed for a flow (a) weakly and (b) strongly deflected. Here, H and D are the crater depth and width, respectively, and L is the distance above the granular bed .

Figure 2 .

 2 Figure 2. Sketch of the device to study the erosion of a granular bed by a water jet: (a) hydraulic circuit to prepare the initial bed by fluidization, and (b) water jet injector dimensions.

Figure 3 .

 3 Figure 3. Experimental results of the crater (a) depth and (b) width normalized by the jet width b as function of the normalized jet-bed distance L/b for different jet Reynolds numbers.

Figure 4 .

 4 Figure 4. Experimental results of the crater (a) depth and (b) width normalized by the reduced jet-bed distance Lλ as a function of E for different jet Reynolds numbers. The solid line corresponds to the linear fit model through the data.

Figure 5 .

 5 Figure 5. Numerical simulations at different stages for the standard two-phase flow model: the solid concentration with fluid field (top), and velocity magnitude of the fluid (bottom).

Figure 6 .

 6 Figure 6. Solid concentration (top) and fluid velocity magnitude (bottom) for the proposed two-phase flow model using a solid-liquid transition model of the sediment phase at different stages.

Figure 7 .

 7 Figure 7. Numerical results of the liquid-solid transition function F at different stages. The fluid velocity field is also plotted as a reference.

Figure 8 .Figure 9 .

 89 Figure 8. Numerical results of the fluid velocity (a) field, and (b) magnitude with streamlines at t = 1.5 s using b = 8 mm, U J = 0.323 m/s, and L = 4.6 cm for different grid resolutions.

  ) and (d)).

Figure 10 .

 10 Figure 10. Velocity field and crater shapes for different jet configurations using b = 4 mm and fine mesh (∆x = ∆z = 0.5 mm).

Figure 11 .

 11 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 models.

Table 2 .

 2 Numerical results of the crater dimensions using different parameters (U J , L, b) and λ = 8.75.

	Mean velocity	Impingement distance	Jet width	Jet-bed distance	Reynolds number	Erosion parameter	Crater width	Crater depth
	U J (m/s)	L (mm)	b (mm)	L/b	Re	E	D (mm)	H (mm)
	0.18	40	4	10	720	2.24	12	3
	0.18	50	4	12.5	720	1.29	15	5
	0.28	50	4	12.5	1120	2.01	32	10
	0.28	60	4	15	1120	1.56	35	6
	0.30	40	4	10	1200	3.74	23	10
	0.30	60	4	15	1200	1.67	40	10
	0.37	46	4	11.5	1480	3.11	40	15
	0.37	80	4	20	1480	1.54	50	10
	0.28	70	6	11.66	1680	2.28	38	15
	0.28	80	6	13.33	1680	1.82	42	13
	0.18	80	8	10	1440	2.24	25	10
	0.18	90	8	11.25	1440	1.58	30	10

Version September 22, 2022 submitted to Water https://www.mdpi.com/journal/water

Funding: This research was partially funded by the Quebec Ministry of Foreign Affairs (MRIF, project AVOCAP-2021-2023), the Mexican Council of Science and Technology (Investigadoras e Investigadores por Mexico, CONACYT). Additional supports by the NSERC-Discovery program (RGPIN-2018-0677) and ComputeCanada (project SINAPSE, No. 3148) are also acknowledged. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.

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.

Author Contributions: Numerical methodology, MUZ, DPVB, WZ, KDN.; experimental study, GG, PG.; writing, DPVB, MUZ. All authors have read and agreed to the published version of the manuscript.

Abbreviations

The following abbreviations are used in this manuscript: