2.4.1. General Setup

The numerical domain dimensions presented in Figure 1 are similar to the ones used by Lee et al. (2016) [12]. For both configurations, the cylinder was placed 15 diameters away from the inlet. The overall domain dimensions were 35 diameters long and 6.1 diameters high. The top boundary condition was a symmetry plane. For the reduced pressure, the outlet boundary condition was a homogeneous Dirichlet condition (*p* − *ρf gy* = 0 Pa). For the outlet velocity, a homogeneous Neumann boundary condition was used for outgoing flows, and a homogeneous Dirichlet boundary condition was used for incoming flows. The inlet was decomposed into two parts. From the bottom to *y* = 1.5 *D*, a wall-type boundary condition was applied. From *y* = 1.5 *D* to the top, a rough wall log law velocity profile was used following the expression:

$$
\mu\_1^f(y) = \frac{\mu\_\*}{\kappa} \ln \left( \frac{30y}{k\_s} \right),
\tag{33}
$$

where *κ* = 0.41 is the von Karman constant, *u*∗ is the friction velocity, and *ks* = 2.5*d* is the Nikuradse roughness length. The different boundary conditions can be found in the test case 2DPipelineScour, publicly available on GitHub (after the paper is accepted). Second order schemes (Gauss linearUpwind) and the default preconditioned biconjugate gradient pressure solver were used for all the simulations presented in this paper.

**Figure 1.** Sketch of the geometry and the boundary conditions used for the computational domain.

For the two configurations, following Chauchat et al. (2017) [17], the turbulent parameter *B* was set to *B* = 1. According to Van Rijn (1984) [27], the value of *σc* depends on the suspension number *ws*/*<sup>u</sup>*∗ with *ws* the particles' fall velocity:

$$\frac{1}{\sigma\_{\varepsilon}} = 1 + 2 \left[ \frac{w\_{\ast}}{\mu\_{\ast}} \right]^2, \ 0.1 < \frac{w\_{\ast}}{\mu\_{\ast}} < 1 \tag{34}$$

Therefore, *σc* is bounded between 1/3 and one. For the present configuration, the sediment transport was intense; the suspension number was small; and following Lee et al. (2016) [12], *σc* was set to one.

The mesh was generated using the OpenFOAM utility snappyHexMesh. Cells were refined in the sediment bed region. Non-refined cells were squares having 3 × 10−3m sides, and refined cells were squares having 7.5 × 10−<sup>4</sup> m sides. The different turbulence models required a specific near-wall resolution. Therefore, cells' refinement close to the cylinder and turbulent boundary conditions depended on the choice of the turbulence model. More details are available in the test case 2DPipelineScour.

#### 2.4.2. Simulations with the *k* − *ε* Turbulence Model

The first cells near the cylinder were 6 × 10−<sup>4</sup> m thick, giving a dimensionless near-wall cell thickness equal to the required *y*<sup>+</sup> = 30 with *y*<sup>+</sup> the dimensionless wall distance defined as *y*<sup>+</sup> = *<sup>u</sup>*<sup>∗</sup>*y*/*<sup>ν</sup>*. A homogeneous Dirichlet boundary condition of 1 × 10−<sup>10</sup> m<sup>2</sup> s<sup>−</sup><sup>2</sup> was applied on the cylinder surface for *k*, and a homogeneous Neumann boundary condition was applied for the rate of dissipation of the TKE *ε*. Similar to Lee et al. (2016) [12], inlet values for turbulent quantities were set following Ferziger (2002) [28]: *k* = 10−4*U*¯ and *ε* = *k*3/2/0.1 h, with h the distance from the bed to the top boundary.

#### 2.4.3. Simulations with the *k* − *ω*2006 and the Modified *k* − *ε* Turbulence Models

For the simulations using the *k* − *ω*2006 and the modified *k* − *ε* turbulence models, cells near the cylinder were 2 × 10−<sup>5</sup> m thick, giving a near-wall cell thickness equal to *y*<sup>+</sup> = 1. Wall functions for smooth walls were applied on the cylinder surface for *k* and *ω*. Inlet values for *k* and *ω* were calculated similarly to Section 2.3.1 following Ferziger (2002) [28] with *ω* = *<sup>ε</sup>*/*Cμk*, except that the dissipation was enhanced at the inlet by two orders of magnitude to reduce the incoming TKE, susceptible to damping the vortex-shedding.
