Background

Total suspended solids (TSS) loads, which is carried out by stormwater runoff during wet weather, is considered as a major pollutant with deleterious effects on receiving water bodies, by transporting adsorbed pollutants, increasing water turbidity, inhibiting plant growth and diversity, blanketing spawning grounds, affecting river biota, and reducing the number of aquatic species [1]. Stormwater detention/retention ponds are one of the most commonly used hydraulic structures in stormwater management both in flooding control by mitigating the peak flow [2,3,4,5], and stormwater runoff quality control by mainly trapping the incoming suspended solids [6,7,8]. TSS trap efficiency of ponds is one of the most informative indicators for evaluating the performance of such facilities and some regulatory agencies have even proposed legislations for sediment removal criteria for stormwater ponds. For example, The City of Calgary has mandated the removal of a minimum of 85% TSS for particle sizes greater than or equal to 50 μm [6]. Despite the significant capacity of removal sediment of the detention ponds, there still exists a gap between the real-world performance and the local legislation or best practice standards for these facilities. Furthermore, the flow and sediment loads entering the stormwater are highly different, resulting in large diversity in removal performance. Birch et al. [9], for instance, reported the range varying from − 12 to 93%, while Gonzalez-Merchan [10] estimated the annual removal efficiency of TSS in Django Reinhardt detention basin located at Lyon in France, varying from 33 to 75%, according to in situ monitoring data during the year of 2004–2010. However, the major influence factors resulting in the significant fluctuation of trap efficiency are still not fully understood [1, 6, 11,12,13,14], such as physical particle characteristics (e.g., sizes and density), pond shape, residence time, vegetation, wind, and turbulence level. As such, it is critical that a comprehensive understanding of various key factors in trap efficiency, thereby, optimizing the design of stormwater ponds to meet the specifications laid out in local legislations. Sedimentation, which includes the process of particulates’ transport, settling to the bottom of a water column, and resuspension due to strong flow, is the dominant mechanism of trap efficiency of ponds; therefore, gaining insight into the sediment transport process under complex conditions is of great importance and is the key difficulty of this topic.

Vertical integrated 2D model based on Saint–Venant equation and 3D models’ coupling with species convection diffusion transport equation has been developed to investigate the sediment transport in stormwater ponds. For example, Gharabaghi et al. [15] and Torres et al. [16] simulated flow and sediment transport of ponds with 2D Rma & Sed 2D and Rubar20, respectively. In their study, the sediment deposition and trap efficiency were modeled; however, variability and the heterogeneous characteristics of sediment, which are essential for removal performance in ponds, have not been considered since both of them only took single particle size into account [6]. To account for the heterogeneous characteristics of sediment in real life, Euler–Lagrange approach has been developed in which the sediment is treated as individual discrete element, actually, over the past decades, and a growing number of researchers have used discrete phase model (DPM) in multiphase model in attending to the sediment transport problems, such as modeling sedimentation in stormwater tank [17,18,19,20,21,22], in storage chamber [23,24,25], in hydrodynamic separator [26] or in Combined Sewer Overflow (CSO) structures [27]. Most of these researches focused on some stable flow patterns [18, 24, 25], which cannot allow for the time-dependent effects. However, as reported by many researches [22, 28,29,30,31,32], unsteady flow is necessary and essential to understand the sediment dynamics. In parallel with the time-dependent condition, boundary condition treatments for DPM model are another important issue to well represent the sedimentation dynamics in which the particle deposited and resuspended are included. However, previous researches on modeling of sediment transport in stormwater pond with DPM have just considered the settling down process, but not the resuspension process due to flow change by time, which is critical for predicting both for trap efficiency and sediment deposition zones [17,18,19,20, 22].

Despite the great progress has been made by the outlined studies in building knowledge for analyzing the sediment transport in stormwater pond, limitations of these numerical models still persist for offering insight into the following: how variability of particle characteristics affects the sediment transport process, how flow pattern alters the predicting of sediment deposition zone, which type of bed boundary conditions and parameters resulting in better modeling performance, what is the role of the flow turbulence in assessing the trap efficiency, and what is suitable particle amount to achieve the target accuracy in terms of computational efficiency. Such detailed understanding of sedimentation process at grain scale is of great importance as it (1) builds a fundamental understanding of sedimentation process, thereby identifying the key impact factors on trap efficiency; (2) provides guidance for practitioners to select appropriate model and parameters for given modeling problems; (3) advances knowledge for optimizing design of stormwater ponds aiming at controlling stormwater runoff quality.

To provide insight into the sedimentation complexity process (transport, settling, resuspension) under time-dependent condition (variable flow pattern and sediment loads), a new model framework based on DPM was proposed, which mainly focuses on the bed load movements at grain scale under unsteady conditions, thereby better understanding of sediment dynamics and more accurate prediction of the sediment on spatial and temporal distributions as well as the trap efficiency in stormwater ponds. The approach was assessed with the experimental data.

Methods

Laboratory experiments

The experiments were designed to investigate the influence of unsteady flow behavior on sediment transport with respect to the settling efficiency and the deposit zones in a shallow tank [22]. The scheme of the experimental setup is shown in Fig. 1. It consists of a rectangular model tank constructed with glass and some adjunctive apparatus. The tank is 3.125 m long, 0.8 m wide and 0.25 m deep. The bottom of the tank has no slope, and a circular pipe inlet (0.1 m diameter) is positioned on the longitudinal centerline at the bottom level of the tank. In the outlet of the model, a weir (weir coefficient µ = 0.62) with an upstream scum board is set to control the water surface level. The inlet tank is required to calm down the inflow in the physical model. The water escapes from the model above the outflow downstream weir. On the bottom of the model, a grid helps to visualize the spatial distribution of the sediment.

Fig. 1
figure 1

Experimental setup

Full size image

To visualize the unsteady flow pattern, and to characterize the shallow tank, the deposition zones and the settling efficiency, both dissolved and particulate tracer tests were carried out. With the dissolved tracer, the visualization of the unsteady flow pattern, especially the position of the inlet jet and the residence time distribution of the shallow tank were determined. For the investigation of the deposition zones, the trap efficiency and the sedimentation processes in the tank, particulate tracer tests were realized.

Sedimentation tests were undertaken with polystyrene particles with a minimum, median and maximum particle size of 300 µm, 500 µm and 800 µm, respectively. The density of the polystyrene particles is 1020 kg/m3. To model the sediment dynamics and the trap efficiency within the tank under unsteady flow, the particles were correspondingly added in the form of impulses through the injection mechanism in 10 steps with a volume of 100 ml at intervals of 5 min. After each experimental test, the volume of particles retained within the tank and the volume of the escaped particles was measured to estimate the trap efficiency. The trap efficiency is defined as the ratio of the sediment retained in the tank to sediment added into the tank. The observations of sediment spatial distribution were made by eye and by photographs at an interval of 5 s during the entire experiment. Thus, the evolution of the sediment transport can be well represented. Three experiments with similar boundary conditions were carried out to examine the stability as a result of the random variations during the experiments.

All measurements and flow simulations for this physical model have previously been reported by Vosswinkel et al. [22]. This study focuses on the development and the assessment of a physically based approach to represent the sediment transport, settling and resuspension processes under unsteady conditions in the tank. The observed data of experiments were then used to validate the numerical modeling.

Numerical methods

Flow modeling

Equations and problem setup

For the incompressible flow in this study, the mass and momentum conservation equations are expressed as Reynolds’ Averaged Navier–Stokes (RANS) form with tensor notation in Cartesian as follows [33]:

The continuity equation

$$\frac{{\partial \left( {\rho \overline{{u_{i} }} } \right)}}{{\partial x_{i} }} = 0.$$
(1)

The momentum equation

$$\frac{{\partial \left( {\rho \overline{{u_{i} }} } \right)}}{\partial t} + \frac{\partial }{{\partial x_{j} }}\left( {\rho \overline{{u_{j} }} \overline{{u_{i} }} } \right) = \frac{{\partial \overline{\tau }_{ij} }}{{\partial x_{j} }} - \frac{\partial P}{{\partial x_{i} }} + \frac{\partial }{{\partial x_{j} }}\left( {{ - }\rho \overline{{u^{\prime}_{j} u^{\prime}_{i} }} } \right) + F_{i} ,$$
(2)

with i and j = 1, 2 and 3, where xi represents the three coordinate axes, \(\overline{{u_{i} }}\) is the time-averaged velocity component in axis i, P means the pressure, ρ is fluid density, \(\overline{{\tau_{ij} }}\) is the mean viscous stress tensor component, \(\overline{{u^{\prime}_{j} u^{\prime}_{i} }}\) = the Reynolds stresses with ‘prime sign’ referring to the time fluctuations, and Fi = the body force (here is the gravity).

The computational meshes play an important role, which pose a significant impact on the solution accuracy, convergence feature and computational time consuming [33, 34]. The independent meshes were tested by Vosswinkel et al. [22] for the same physical model in the previous flow modeling. Accordingly, the same mesh was used in this study. A preliminary study has showed that the k-ε RNG turbulence model can well reproduce the unsteady flow patterns, particularly for the periodical oscillations flow, comparing to observations [22]. Hence, the k–ε RNG model was also beneficial to unsteady sediment transport modeling in this study.

Sediment transport modeling with DPM

Basic equations’ description

The Lagrangian DPM is derived from force balances based on Newton’s law describing particle trajectory. The force balance in the x-direction (in Cartesian coordinates) can be written as [35]

$$\frac{{{\text{d}}u_{p} }}{{{\text{d}}t}} = F_{D} \left( {u - u_{p} } \right) + \frac{{g_{x} \left( {\rho_{p} - \rho } \right)}}{{\rho_{p} }} + F_{x}$$
(3)

where up is the x-velocity component of the particle, u the instantaneous x-velocity component of the fluid, gx the gravity acceleration along the x-axis, ρp the density of the particle, ρ the fluid density, and Fx the additional forces, such as the Saffman’s Lift Force due to the shear effect or other forces dependent on the case [35]. The left term in Eq. (3) corresponds to the x-direction acceleration of the particle per unit particle mass. The first term on the right-hand side corresponds to the mass force of drag per unit particle mass. The coefficient FD is expressed by Eq. (4).

$$F_{D} = \frac{18\mu }{{\rho_{p} d_{{}}^{2} }}\frac{{C_{D} {\text{Re}}_{p} }}{24},$$
(4)

where μ is the dynamic viscosity of the fluid, d the diameter of the particle, CD drag coefficient, and Rep the particle Reynolds number. The second term comprises forces due to gravity per unit particle mass (competition between weight and buoyancy of the particle). The third term corresponds to additional forces per unit particle mass, including the added mass force and the force due to pressure gradient. Rep is a dimensionless number characterizing the relative velocity between fluid and particle, defined in Eq. (5):

$$\text{Re}_{p} \equiv \frac{{\rho d_{p} \left| {u - u_{p} } \right|}}{\mu }$$
(5)

where \(u\) and \(u_{p}\) are flow and particle velocity, respectively. The drag coefficient CD depends on the flow regime. A drag coefficient expression developed by Morsi and Alexander [36] for spherical particles, which provides the most complete range of Rep, was widely used [18, 25, 26]. This expression is written as

$$C_{D} = \alpha_{1} + \frac{{\alpha_{2} }}{{\text{Re}_{p} }} + \frac{{\alpha_{3} }}{{\text{Re}_{p}^{2} }}.$$
(6)

The empirical constants α1, α2 and α3 depend on the particle Reynolds number [36].

Turbulent dispersion effect

Some researchers [37] assumed that particle trajectory can be calculated with the mean flow without considering the influence of flow turbulence. Others [18, 19, 22, 25, 38,39,40,41,42] argued that the dispersion of small particles is severely affected by flow turbulence.

Dufresne et al. [18] and Vosswinkel et al. [22] used the stochastic approach (discrete random walk model—DRWM) to represent the turbulent dispersion effect and tested the parameter concerned, the integral time scale constant CL. A large CL value means that the particle trajectory is strongly affected by the turbulence. These authors found that the simulated trap efficiency and sediment deposit zone were sensitive to this parameter. In this study, turbulent dispersion effect was modeled by DRWM. The DRWM contributes to the effect of instantaneous turbulent velocity fluctuations on particle trajectories through stochastic methods by eddies. Technically, each eddy is characterized by a Gaussian distributed random velocity fluctuation u’, v’, and w’. Velocity fluctuations are expressed in Eq. (7), where ζ is a normal distributed random number while k is the turbulent kinetic energy.

$$u^{\prime} = v^{\prime} = w^{\prime} = \zeta \sqrt {\frac{2k}{3}}$$
(7)

In the DRWM approach, the fluctuating velocity components are discrete piecewise constant functions of time. The time scale constant CL within the DRWM is not well known and was used as a calibrated parameter in the previous studies [22]. For more details on DRWM, one can refer to Quadrio and Luchini [43], Dufresne et al. [18] and Ansys [35].

Sediment particle size distribution

The Rosin–Rammler distribution function was employed to represent the non-uniform sediment sizes. The whole range of particle sizes is divided into a number of discrete intervals. Each interval is represented by a mean diameter and tracked for the particle trajectory calculation. The Rosin–Rammler distribution (developed by Rosin and Rammler [44]) is an empirical distribution to describe particle sizes distribution, which can be expressed as

$$Y_{d} = 100\exp \left[ { - \left( {{\raise0.7ex\hbox{$d$} \!\mathord{\left/ {\vphantom {d {\overline{d} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\overline{d} }$}}} \right)^{n} } \right]$$
(8)

where Yd is cumulative mass in % retained on size d,\(\overline{d}\) the size parameter, and n a distribution parameter. The parameter of Rosin–Rammler model is determined by the measured data, and more details can be found in Ansys [35]. The setting details of Rosin–Rammler model are listed in Table 1.

Table 1 Basic configurations for DPM model excluding bed boundary conditions
Full size table
Treatment of the interaction between particles and the bed under unsteady condition

With DPM modeling, two kinds of boundary conditions are often used in earlier investigation to deal with the particle trajectory fate when hitting the bed: trap and reflect [18, 25].

The trap boundary condition means that particles are trapped and end its tracking calculation when they arrive at the physical boundary, such as the bed of the tank. This boundary condition excludes the possibility of particle resuspension after hitting the bed and therefore often overestimates sediment trap efficiency in a storage tank [23]. The reflect boundary condition means that the particle rebounds off the boundary via an elastic collision, which excludes the possibility of the particle settling out when it arrives the bed. Therefore, when acting as the boundary condition at the bed of tank, the reflect boundary condition often underestimates sediment trap efficiency in sediment transport modeling in a storage tank [24, 25].

Adamsson et al. [25] developed a boundary condition for the bed of the tank by combining these two basic boundary conditions (trap and reflect) together with critical bed shear stress (BSS). If the local bed shear stress is less than the critical value, the particle settles out as with trap; otherwise, the particle rebounds off the bed as with reflect. Likewise, Dufresne et al. [18] used bed turbulent kinetic energy (BTKE) as the threshold in the same way to treat particle trajectory fate when it hits the bed of tank. However, Adamsson et al. [25] pointed out that this boundary treatment was only valid for steady flow condition and the prediction was sensitive to the chosen threshold. It cannot represent the dynamic morphology evolution of sediment such as resuspension process under unsteady condition. That is why the model prediction often overestimates the trap efficiency even for large inflow rate during experimental tests [18, 24].

To accurately model the dynamic sedimentation and resuspension phenomena, an approach aiming at representing the dynamic settling of suspended sediment and resuspension of deposit under unsteady conditions based on DPM is presented. The approach is shown in the diagram in Fig. 2. In addition, to accurately reproduce the flow and sediment transport behaviors, the liquid and solid phases are integrated in the two-way coupling which takes the interaction between flow and sediment into account for their important influence as reported by Kantoush et al. [30] and Stovin and Saul [45].

Fig. 2
figure 2

Diagram of dynamic sedimentation and bed load resuspension under unsteady flow conditions

Full size image

The particle’s state which varies from being suspended load, bed load and deposit are controlled by the near bed flow conditions. As shown in Fig. 2, at the grain scale, two sediment states are clearly distinguished at the moment t: it is moving (suspended load or bed load) or settled out as deposit. Here just those moving particles which are going to hit the bed and the deposit need to be checked for state change. The state of a particle at the moment t and the next moment t + Δt may change (deposit or moving, corresponding to the settling process and the resuspended process, respectively) due to the flow variability. The treatment approach for these particles shown in Fig. 2 is described as follows.

Process A: particle settling For those particles in motion at the moment t, there are two possible states at the next moment t + Δt, either settling out as deposit (condition A1 in Fig. 2) or continuing to move as suspended load or bed load after contacting with the bed (condition A2 in Fig. 2). The bed shear stress (BSS) is used to determine the state of the particle after contacting with the bed. The Shields curve [46] was used to calculate a BSS threshold according to the particle characteristics. Equations (9) and (10) show a set of fitted formula of the Shields curve, termed as VBSS (varying bed shear stress) hereafter. The critical bed shear stress threshold then can be estimated by Eq. (11). If the local BSS value is lower than the threshold, the particle settles down as deposit, or else it continues to move after contacting with the bed.

$$\tau_{*c} = 0.22\beta + 0.06 \times 10^{ - 7.7\beta }$$
(9)
$$\beta = \left[ {\frac{\rho }{\mu }\sqrt {\frac{{\rho_{\text{p}} { - }\rho }}{\rho }gd^{3} } } \right]^{ - 0.6}$$
(10)
$$\tau_{c} = \tau_{*c} \left( {\gamma_{s} - \gamma } \right)d$$
(11)

where τ*c is the dimensionless shear stress, as is particle specific weight, γ is fluid specific weight, μ is dynamic viscosity of fluid, g is acceleration of gravity, d is particle diameter, pp. and ρ are particle and fluid density, respectively.

Condition A1: If the particle settles down from motion, the velocity and the forces associated with particle movement should be adjusted so it actually remains as a deposit. In this way, the deposit may be resuspended if surrounding conditions change.

Condition A2: If the particle still keeps on moving after coming into contact with the bed, a variant reflect boundary condition is preferable. The particle motion decomposition is shown in Fig. 3. The striking particle velocity is decomposed into normal and tangential components with respect to the collision surface, vin and unit, respectively. The rebound particle velocity components are under control of the normal coefficient e and tangential coefficient f, respectively [47, 48]. Technically, the normal coefficient defines the amount of momentum in the direction normal to the boundary after collision. By contrast, the tangential coefficient identifies the amount of momentum in the tangential direction retained by the particle after collision. A normal and tangential coefficient of restitution up to 1 implies that the particle retains all of its normal and tangential momentum (no energy lost) after the rebound. On the contrary, a normal or tangential coefficient of restitution equal to 0 implies that the particle retains none of its normal or tangential momentum after the rebound (particle energy lost entirely). For the condition of e equal to 0, the sliding or rolling behavior is represented. An empirical model involved in the restitution coefficient and tangential friction resistance coefficient has been put forward by Niño and García [47] with Eq. (14). The particle moving direction is assumed to keep the same particle longitudinal direction before the collision.

Fig. 3
figure 3

Diagram of the rebound boundary condition

Full size image
$$e \equiv \left| {\frac{{v_{out} }}{{v_{in} }}} \right|$$
(12)
$$f \equiv \frac{{u_{out} }}{{u_{in} }}$$
(13)
$$\left\{\begin{array}{*{20}l} {e = 0.84 - 4.84\tau_{*}} \hfill \\ {f = 0.73} \hfill \\ \end{array} \right.$$
(14)

where τ* is the dimensionless shear stress, e and f indicate the normal and tangential coefficient, respectively.

Process B: particle resuspension The deposit may be resuspended as bed load or suspended load at the moment t + Δt with flow conditions’ change. The local bed shear stress is compared to the threshold of shear stress calculated by Eqs. (9)–(11) according to the particle characteristics. If the local shear stress is larger than the threshold, the particle is entrained (condition B2 in Fig. 2); or else, the particle remains in the same location as deposit.

Condition B2: If the deposit is entrained, an initial velocity should be specified. In this study, the formulas developed by Hu and Hui [49] described by Eqs. (15)–(16) are adopted for the estimation of the initial velocity:

$$u_{i} = \left\{ {\begin{array}{*{20}l} { \left( {12.3 - 3.7\log \tau_{*} } \right)u_{*},} \hfill & {if\,\,\tau_{*} < 1.2} \hfill \\ {12.1u_{*} ,} \hfill & {if\,\,\tau_{*} < 1.2} \hfill \\ \end{array} } \right.$$
(15)
$$v_{i} = \left\{ {\begin{array}{*{20}l} { \left( {3.2 - 4.5\log \tau_{*} } \right)u_{*} ,} \hfill & {if\,\,\tau_{*} < 1.2} \hfill \\ {3.1u_{*} ,} \hfill & {if\,\,\tau_{*} > 1.2} \hfill \\ \end{array} } \right.$$
(16)

where ui and vi are the tangential and normal velocity components for the particle at the hitting surface, respectively (m/s), u* = (τ0/ρ)0.5 is the shear velocity (m/s), and τ* is the dimensionless shear stress.

Regarding the particle velocity direction, it is reasonable to assume that the resuspended particle follows the same longitudinal direction as the fluid flow, since the particle is entrained by the near bed flow. The developed approach for dynamic settling and resuspension is implemented with an additional module with the user-defined functions (UDF). The UDF is coded in C language, and then is compiled in Fluent and hooked into Ansys Fluent version 14 as a bed boundary condition of the tank for unsteady DPM modeling under unsteady flow conditions. The developed UDF code for the developed approach is provided in the supplemental file, and the detailed compile and setup procedure can refer to the official document [50].

Model setup

The flow and DPM modeling settings are given in Tables 1, 2 and 3. The basic DPM configurations excluding the bed boundary conditions are listed in Table 1. Table 2 shows the main settings of fluid flow modeling. The developed approach is specified as boundary condition in the bed of tank via UDF for DPM to represent the dynamic sediment transport behaviors near the bed (settling, sliding, saltating, and resuspension). To find the optimal parameters for the developed approach, a series of tests with different settings for the approach were designed systematically. The different parameter setups are shown in Table 3.

Table 2 Configurations of flow modeling
Full size table
Table 3 DPM setup and different bed boundary conditions for the developed approach
Full size table

Trial and error tests were undertaken with different quantities of injected particles for sediment transport modeling to ensure the statistical characteristics with an appropriate quantity of particles for reducing the computational resources’ requirement. The test cases DPM 1 to DPM 5, shown in Table 3, were designed for this purpose. The cases DPM 6–DPM 9 and DPM 11 shown in Table 3 were designed to assess the influence of the turbulence dispersion by means of DRWM. Cases DPM 10 and DPM 11 were designed to evaluate the influence of the resistance coefficient on the sediment dynamics behavior. Case DPM 12 was designed to evaluate the influence of the resuspension model in the developed approach. The validation of CFD models was subjected to comparisons of sediment transport predictions against experimental data. Validation involved the qualitative comparison of experimental sediment in temporal and spatial distribution with numerical prediction and the quantitative comparison of measured and simulated trap efficiency.

Although many test cases were carried out, it should be noted that some assumptions are involved in the CFD models. They are summarized as follows: (1) the particles were considered as a sphere; (2) the particles’ density is unique; (3) the particle size distribution is much like a Rosin–Rammler distribution; (4) the accumulation of sediment on the bed is thin and spare enough so that the bed reflect condition is reasonable; (5) the particle in the model represents a parcel of the real particles which have the similar physical properties and total number of injected particles is not the exact number of the real condition, but with the same total injected mass or concentration in the experiment.

Results and discussions

The visual comparison would give a sight into the potential of the developed approach in representing sediment transport under unsteady-state conditions. Then, the further advantage of this approach is analyzed. Finally, the parameters or models involved in the approach are primarily discussed for best fitting with experimental data.

Flow patterns

The visualization with impulse-like fluorescent tracer (one can refer to Vosswinkel et al. [22]) and particles indicated that the flow pattern developed in the tank was highly unstable although the inflow rate was steady for the experiments, in particular, the recirculation of the inlet jet and the sweeping of the jet are indicators for the unsteady flow pattern in the tank, and the similar results were also reported by others [18, 32, 33, 51]. These results highlighted the time-dependent effect on the sediment transport processes, which is different from the results of the experiments reported by Adamsson et al. [25] and Dufresne et al. [18] for which the flow patterns were steady. Observation shows that the particles were driven by flows and the deposit particles were entrained due to the unstable flow conditions.

According to the experimental observation, the inflow jet first shifted to the left side of the tank to form a counter-clockwise circulation, sweeping from the half length of the tank to the upstream left corner of the tank. Then, the circulation detached from the left side wall and shifted to the right-side wall of the tank, forming a clockwise circulation and sweeping from the half length of the tank to the upstream right corner of the tank. The circulation went back to the left side again in a short period of time. The same experiments were repeated three times and similar processes were observed. The inflow jet oscillation from side to side is subjected to a quasi-periodic behavior. These unstable flow patterns also highlighted the necessity and importance of the time-dependent effect on the sedimentation transport processes, which is hard to represent with the steady-state model developed by the previous researches [18, 24, 25].

Comparison of experimental observation and sediment transport modeling

Due to the limitations of measured devices, it was hard to obtain precise sediment spatial distribution data to verify the modeling results quantitatively. Here, we employed the successive (every 5 s) photograph data to first assess the numerical modeling results by comparing the experimental sediment transport morphology with the simulated results. The simulated results (the sediment morphology) were extracted at the similar viewpoint and direction and compared with the observed data. Figures 4 and 5 display the experimental and simulated evolution of sediment transport processes at the beginning stage and right circulation stage. Both of them demonstrated the reasonable agreement between the simulated results and the experimental observations on the spatial distribution of particles and preferential profile, however highlighted the different nature of the sediment dynamics. Whereas the first one demonstrated the suspended particles’ advection process, the second one indicated the complex bed load movement, which means that it involved particles’ transport, settling down and resuspended processes. In particular, Fig. 5 highlights the most active region of bed load motions (as shown with dotted line both in observation and simulated results). It can be clearly observed that the simulation reproduces the flow circulation evolution pattern (with the clockwise arrow in the observation and velocity arrow in simulated results). Like the flow pattern, the pattern of the bed load was evaluated similarly including suspension settling down as deposit and deposit entrained as suspension. More details will be discussed in “Analysis of the simulated sediment transport processes” section.

Fig. 4
figure 4

Comparison between the experimental and simulated evolution of sediment transport: (above) simulated morphology of sediment in volume concentration at 10 s and 15 s, respectively (DPM 10); the dotted lines indicate the sediment profile at the similar view to the experimental observation; and (below) observed morphology of sediment after particle injection for 10 s and 15 s, respectively; the dotted lines indicate the sediment profile from the simulated results at the same time for comparison

Full size image
Fig. 5
figure 5

(Above) observed flow and sediment patterns evolution at the moment of 90–105–120 s from left to right, respectively, with dotted line and curved arrow highlight the inflow jet clockwise circulation (near the right-side upstream corner) [22]; and (below) simulated flow and sediment patterns evolution at the moment of 90–105–120 s from left to right, respectively, with velocity arrow highlight the inflow jet clockwise circulation (it should be noted that the velocity arrows indicate the largest flow velocity in the whole domain, the lower ones in other zones were filtered) and solid line indicated the preferential shape of bed load distribution from the observation for comparison

Full size image

Moreover, the preferential sediment distribution in the downstream part and in the left and right-hand upstream corners of the tank was also predicted as shown in Fig. 5. Comparison between the simulated and experimental results, it seems that the observed active sediment zones (with dotted line) were closer to the upstream wall than the simulated active zones (with solid line); this is considered to be resulted from the flow field simulation, but not the developed sediment transport model scheme. Due to the measurement limitation, it is hard to obtain more precise measured data to assess the verified simulation. Nevertheless, it can be concluded that the developed approach is able to reproduce the sediment transport processes including advection transport, settling and resuspension processes under unsteady conditions with some model assumptions.

Analysis of the simulated sediment transport processes

Sediment transport micro-mechanism

Different from the previous sediment transport modeling in tank under steady flow condition [18, 24, 25], this study focused on the sediment dynamics under unsteady conditions. Figure 6a displays the numerically simulated deposit spatial distribution after 105 s since the first injection of sediment. Figure 6b demonstrates the deposit and the near bed moving sediments (bed load), which is within a thin layer less than 0.002 m from the bed. Comparison between Fig. 6a and 6b indicated the dynamic sediment transport characteristics. It illustrated that the deposits around the clockwise recirculation were eroded due to the strong flow vortex (large flow velocities around the circulation currents) and were moving with some relatively large velocities. The trend in bed load movement is clearly represented, indicating a reasonable level of agreement with the experimental photographic observations. Similarly, the suspended particles can be represented in the same way. Form Fig. 6b, it can be also observed that the particles in the downstream part of tank were quite inactive because of the low flow velocities. This can explain why most particles were accumulated in that region. This modeling is not only able to the prediction of sediment distribution, but also able to provide information to analyze the sediment dynamics. This means that the verified modeling data would be a great complementary database for sediment dynamics analysis, which were very hard and expensive or even impossible by laboratory experiments and observation in situ. This is a marked advantage compared to the previous related researches work [25], which just establishes the final simulated results such as the sediment zone and trap efficiency for certain steady flow patterns. Based on the developed unity of methods, the simulated results are able to provide detailed information of the sediment to gain an insight into the sediment transport mechanisms at grain scale and thus to better understand the sediment transport dynamics. Of course, it should be noted that the models must be verified carefully as accurate as possible before adopting the simulated data for analysis since some assumptions have been made for model construction. As mentioned in “Methods” section, in DPM model, each spherical particle represents a parcel of particles with similar properties; therefore, the simulation demonstrates a comprehensive and representative result, but not an exact real-world condition. In addition, interaction between particle and particle is ignored in DPM model, which is reasonable in spare particle transport condition. However, this DPM model may lead to unacceptable error in multi-phases systems with dense particles’ condition.

Fig. 6
figure 6

Comparison of sediment spatial distribution out of simulation and experiments. a Simulated deposits in the tank (at the moment of 100 s, the legend indicates the different particle size in meter); b simulated deposits + moving bed load with a thin layer (less than 0.002 m), with the arrow indicating the velocities of the particles and the dotted line indicating the observed preferential sediment zone; and c overhead view of experimental sediment spatial distribution after 100 s of first injection of sediments, with the dotted line indicating the observed preferential sediment zone for the purpose of comparison in the most active zone and the curved arrow indicating the vortex direction

Full size image

Trap efficiency

In terms of experimental trap efficiency, three tests were performed in the physical model. The overall trap efficiencies were 83%, 85% and 87%, respectively, with an average efficiency of 85%. These results were higher than those reported by Stovin and Saul [23, 24] and Dufresne et al. [18], probably due to the different outlet configuration. The previous experiments’ outlet configurations were mostly pipe or narrow channel with inverse altitude near the bed the tank, whereas the present experimental outlet was set as overflow weir with a scum board, which would stop the particles escaping from the low bed level. These trap efficiency prediction models based on the measured data were not suitable for the present case any longer. All of them underestimated the trap efficiency for the present experiments. In fact, Dufresne et al. [51] pointed out that it needs to use different regression laws depending on the flow pattern. Indeed, the relationship laws appear to be case-specific since they were based on some limited test conditions (geometry of tanks, sediment characteristics and flow patterns).

Compared to the conventional regression models, the developed approach is more robust and universal because it is based on the fundamental principle. Moreover, it has accounted for the influence factors on geometric properties of tank, variability of hydraulic and characteristics of sediments. Therefore, it has the capability of predicting the global efficiency without limitations like the regression models. Table 4 presents the simulated trap efficiencies with different parameter settings in the developed approach. All the cases except for case DPM 9 and 12 predicted reasonable trap efficiencies. Case DPM 9 largely overestimates the trap efficiency due to no consideration of turbulent dispersion, while DPM 12 owing to no consideration of resuspension process. Besides, it should also be noted that the trap efficiency is sensitive to the total number of injected particles due to the particle parcel assumption according to the DPM cases 1, 2, 3, 4 and 5. Overall, despite its good capability of high accurate trap efficiency prediction, model parameters should be carefully selected for good trap efficiency prediction since it seems to be affected by many factors. Further discussions are shown in “Influence of resuspension model under unsteady condition” section.

Table 4 Simulated trap efficiencies of different case setups
Full size table

Sensitivity of parameters in the developed approach

Dispersion effect of turbulence on sediment transport

As mentioned previously in “Sediment transport modeling with DPM” section, arguments are existent for the influence of turbulence on the particle motion. This study tested one case without DRWM and several cases with DRWM in different parameters for calibration. Figure 7a, b shows the different sediment spatial distributions with or without stochastic modeling (DRWM), respectively. Compared to the experimental observation shown in Fig. 7c, the simulated results with DRWM show a higher level of agreement than that without DRWM (Fig. 7b). For example, there was a missing zone in the simulated results without DRWM (shown as ellipse solid line in Fig. 7b) as well as the uncovered zones in the two upstream corners of the tank. Whereas the simulated results with DRWM showed a more reasonable sediment spatial distribution across the tank (Fig. 7a). For the trap efficiency prediction shown in Fig. 8, the simulated result without DRWM (DPM 9) obviously overestimates the trap efficiency, with 96.65% against the averaged measured efficiency around 85%. Therefore, both the analyses of the sediment distribution and the trap efficiency may suggest that sediment transport using DPM in shallow tanks should take turbulent dispersion effect into consideration to represent sediment transport behaviors more accurately.

Fig. 7
figure 7

Comparison of numerically simulated sediment deposition pattern distribution a with stochastic modeling (DRWM, CL = 0.05), b without stochastic modeling (no DRWM), c experimental observation

Full size image
Fig. 8
figure 8

Comparison of simulated and experimental trap efficiencies (the bracket indicates the corresponding parameter settings)

Full size image

In fact, Dufresne et al. [18] used DPM with DRWM to predict deposit zone in a sewer detention tank under steady condition, finding that a larger CL (CL = 0.5–1) can give a better result than that with the default value (CL = 0.15). Whereas Vosswinkel et al. [22] indicated that a smaller CL (CL = 0.12) for DRWM showed a better agreement comparing to the measured residence time distribution profile. It seems that the CL was used more likely as a calibration parameter in those investigations. Up to day, there is no universal CL value or law for all cases in literature. In the present study, out of the five tested CL values (0.05, 0.1, 0.15, 0.25, and 0.5) with DRWM (settings shown in Table 3) under unsteady conditions, Fig. 8 shows that the smaller CL values 0.05 and 0.1 predict the closest trap efficiencies (DPM6: 87.85%, DPM 10: 86.77% and DPM 11: 83.94%) compared to the measured data (83–87%). However, this CL range found in the present research was different from the previous researches [18]. It was difficult to tell which one is correct because it is not well known and may depend on the study case configuration (e.g., total injected particles, steady or unsteady simulation mode). As shown in Fig. 8, trap efficiency prediction changes when the total particle number varies even with the same CL value (DPM 1 to DPM 5). In this study, tests were performed under unsteady conditions, which is closer to realistic condition. Thus, it is considered that the range of CL = 0.05–0.1 is the most appropriate interval for DPM to represent turbulent dispersion.

Influence of resuspension model under unsteady condition

Figure 9a displays the simulated final sediment distribution on the bed of the tank based on the developed approach. Compared to the experimental observation (see Fig. 9b), it shows a better agreement of deposit distribution prediction than the two different, earlier deposit predictions without consideration of deposit resuspension (see Fig. 9c, d). Figure 9c shows the worst prediction because the prediction was based on a steady flow simulation, which indicates that steady flow mode cannot represent the observed unsteady flow behavior. Figure 9d illustrates the superposition of deposit predictions with steady DPM modeling based on three instantaneous flow fields of the unsteady simulation, taken as snapshots of characteristic flow field situations at 650, 750 and 850 s of simulation time, respectively. This analysis technique may be able to predict the deposit zones at different moments, but it cannot account for the time-dependent effect of deposit rearrangement and resuspension, like the presented refined approach does (Fig. 9a). Finally, Fig. 9e displays the simulated deposit distribution based on the developed approach, but without consideration of the resuspension effect. It can be observed that a large amount of sediment settled down in the upstream part and the middle part of tank. In fact, there was much less sediment observed from the experiments in the same zone (see Fig. 9b).

Fig. 9
figure 9

Comparison of a simulated spatial deposit distribution at the end, with consideration of deposit resuspension under unsteady conditions (here proposed refined approach; case DPM 10 at 300 s), b observed final deposit distribution (Vosswinkel et al. [22]), c simulated deposit distribution based on a steady flow simulation, without consideration of resuspension (Vosswinkel et al. [22]), d superposition of three simulated deposit distributions, without consideration of resuspension, in the temporal sequence of 650 to 850 s (Vosswinkel et al. [22]), e simulated final deposit distribution without consideration of deposit resuspension (DPM 12 at 300 s)

Full size image

In addition, without consideration of resuspension the modeling overpredicted the trap efficiency (92.77% of simulation of DPM case 12 against 85% of experiment observation shown in Fig. 8), while a linear resuspension formula gave the best trap efficiency prediction (86.77% and 83.94% of simulation of DPM case 10 and 11 against 85% of experiment observation shown in Fig. 8). This reveals that resuspension effect plays an important role in sediment transport processes, and thus greatly affect the simulated results on sediment distribution and trap efficiency under unsteady conditions. The results show the advantage of the developed approach and suggest that the resuspension model should be considered unconditionally in sediment transport modeling.

Overall, the newly developed unity of numerical modeling approach demonstrates a quite attractive capability of representing the particle transport, settling and resuspension behaviors under unsteady conditions. However, the model prediction accuracy is affected by the involved assumptions and model parameters. An optimal set of parameters was determined with a series of laboratory-scale experiments. Application of the developed approach in engineering still needs further parameter validation since the sedimentation phenomenon in real-world big tanks is far more complex than the laboratory-scale experimental tests. Indeed, it is by no means easy to draw a definitive and universal conclusion as the combination of these parameters is highly complicated and their interaction with each other is still not quite clear.

Conclusions

In this study, a novel modeling framework with a physically based DPM was presented, aiming at better representing the sediment dynamics under unsteady-state conditions. With a range of different parameters involved in the developed approach, a series of different case tests were designed, conducted and analyzed to assess its validity. The simulated results were compared with experimental data qualitatively and quantitatively. The following conclusions were summarized through simulations and observations:

  1. 1.

    The proposed method can well reproduce the evolution of sediment transport, settling and resuspension behaviors with more detailed information than the experimental case at grain scale, it may imply that the decomposition of bed load motion at grain scale with the proposed modeling framework is a valid way to represent the sediment transport phenomena;

  2. 2.

    The proposed method can accurately predict the trap efficiency and temporal–spatial sediment distribution with calibrated parameters;

  3. 3.

    Time-dependent effect is critical and necessary to model sediment transport process, especially for accuracy of deposition zones prediction;

  4. 4.

    Both the turbulent effect with DRWM and resuspension effect have major impact on trap efficiency prediction performance.

The developed model could be a helpful tool for further understanding the sediment transport phenomena at grain scale in shallow tanks since it can provide more detailed information than the experimental tests’ measurement. It is also expected that this model is able to extend easily to more complex phenomena related to suspend particles since other physical models (such chemical or biological models) can be added to these individual particle objects. However, it should also be noted that the developed model is currently not suitable for the multiphase system with complex particle shape and dense particles’ condition due to the spherical particle parcel assumption and non-consideration of particle–particle interaction. Further research would be performed in the next step for more complex system.