### 1. Introduction

### 2. Material and Methods

### 2.1. Background and Mathematical Models

### 2.1.1. Computational fluid dynamics

*k – ɛ*turbulence model, were solved using FLOW-3D software (Flow Science Inc., Santa Fe, NM, USA) to simulate turbulent fluid motion within a retention pond. In the CFD-DPBE model, particles/flocs are assumed to travel via fluid motion and to aggregate or disintegrate due to impact and shear forces or effects [14, 17, 21].

##### (2)

$$\frac{\partial \langle {U}_{i}\rangle}{\partial t}+\langle {U}_{j}\rangle \frac{\partial \langle {U}_{i}\rangle}{\partial {x}_{j}}+\frac{\partial \langle {u}_{j}{u}_{i}\rangle}{\partial {x}_{j}}=-\frac{1}{\rho}\frac{\partial \langle p\rangle}{\partial {x}_{i}}+\nu {\nabla}^{2}\langle {U}_{i}\rangle $$##### (3)

$$\langle {u}_{i}{u}_{j}\rangle =\frac{2}{3}k{\delta}_{ij}-{v}_{T}\left(\frac{\partial \langle {U}_{i}\rangle}{\partial {x}_{j}}+\frac{\partial \langle {U}_{i}\rangle}{\partial {x}_{i}}\right)$$*i*and

*j*are indices,

*x*

*represents coordinate directions (*

_{i}*i =*1, 2, 3 for

*x*,

*y*,

*z*directions, respectively), 〈

*U*

*〉 is the time-averaged velocity component (m/s),*

_{i}*t*represents time (s),

*ρ*is the fluid density (kg/m

^{3}),

*p*is the piezometric pressure (kg/m/s

^{2}), and

*ν*is the kinematic viscosity of the fluid (m

^{2}/s). A symmetric second-order tensor 〈

*u*

_{i}*u*

*〉 representing Reynolds normal or shear stresses (m*

_{j}^{2}/s

^{2}) is modeled with Eq. (3), and

*ν*

*, the turbulent viscosity (N·s/m*

_{T}^{2}), is specified by Eq. (4). In Eqs. (3) and (4),

*δ*

*is Kroneker’s delta,*

_{ij}*C*

*is 0.09, a model constant,*

_{μ}*k*represents turbulent kinetic energy (m

^{2}/s

^{2}), and

*ɛ*is the turbulent energy dissipation rate (m

^{2}/s

^{3}). These two energy terms are obtained from Eqs. (5) and (6), which together represent the

*k–ɛ*turbulence model. They are solved simultaneously with the RANS equations. The velocity gradient ( $G=\sqrt{\varepsilon /v}$, /s), which is obtained from the two-equation

*k – ɛ*turbulence model, causes particle/floc aggregation or break-up kinetics in the DPBE and thus serves as a coupling term between the turbulent flow field (CDF problem) and the DPBE [14]. In Eqs. (5) and (6), model fitting constants have been found to be

*σ*

*= 1.0,*

_{k}*σ*

*= 1.3,*

_{ɛ}*C*

*= 1.44, and*

_{ɛ1}*C*

*= 1.92 [17].*

_{ɛ2}**P**, specified by Eq. (7), is the turbulent kinetic energy production term (m

^{2}/s

^{3}).

##### (5)

$$\frac{\partial k}{\partial t}+\langle U\rangle \xb7\nabla k=\nabla \xb7\left(\frac{{v}_{T}}{{\sigma}_{k}}\right)+P-\varepsilon $$### 2.1.2 Multi-dimensional discretized population balance equations

##### (8)

$$\begin{array}{ll}\left[\frac{\partial {n}_{i}}{\partial t}\right]\hfill & (\text{I})\hfill \\ +\left[\frac{\partial}{\partial x}(\langle {U}_{x}\rangle {n}_{i})+\frac{\partial}{\partial y}(\langle {U}_{y}\rangle {n}_{i})+\frac{\partial}{\partial z}(\langle {U}_{z}\rangle {n}_{i})\right]\hfill & (\text{II})\hfill \\ -\left[\frac{\partial}{\partial x}\left({C}_{\mu}\frac{{k}^{2}}{\varepsilon}\frac{\partial {n}_{i}}{\partial x}\right)+\frac{\partial}{\partial y}\left({C}_{\mu}\frac{{k}^{2}}{\varepsilon}\frac{\partial {n}_{i}}{\partial y}\right)+\frac{\partial}{\partial z}\left({C}_{\mu}\frac{{k}^{2}}{\varepsilon}\frac{\partial {n}_{i}}{\partial z}\right)\right]\hfill & (\text{III})\hfill \\ +{(agg/break)}_{i}-{u}_{gi}\frac{\partial {n}_{i}}{\partial z}\hfill & (\text{IV})\hfill \end{array}$$*n*

*=*

_{i}*n*(

*x, y, z, D*

_{i}*, t*) is number concentration of flocs (/m

^{3}) of linear class size

*D*

*(m) (*

_{i}*i=1, 2, …i*

*;*

_{max}*D*

*≤*

_{1}*D*

*≤*

_{i}*D*

*; for all*

_{max}*D*

*,*

_{i}*n*

*is called the population density function);*

_{i}*x, y, z,*and

*t*are position and time; 〈

*U*

_{x}_{〉}, 〈

*U*

*〉, and 〈*

_{y}*U*

_{z}_{〉}mean fluid velocity components in the

*x, y,*and

*z*directions (m/s);

*ρ*is fluid density (kg/m

^{3});

*k*=

*k(x, y, z, t)*is turbulent kinetic energy (m

^{2}/s

^{2});

*ɛ = ɛ*(

*x, y, z, t*) is turbulent energy dissipation rate (m

^{2}/s

^{3});

*C*

*= 0.09 is standard value of a CFD model constant (see Eq. (4)); and*

_{μ}*u*

*is settlement velocity (m/s) of the*

_{gi}*i*th floc class due to gravity. On the left-hand side of Eq. (8), the respective terms in brackets represent the storage change (I), particle/floc mean advection (II), and the turbulent dispersion of the particle/floc (III), while on the right-hand side, the source/sink terms (IV) represent the net effects of aggregation, break-up, and settling due to gravity [14]. The coefficients or functions in these later terms are largely empirical and must be determined by experiment. The quantities depending on the turbulent fluid variables (〈

*U*

*〉, 〈*

_{x}*U*

*〉, 〈*

_{y}*U*

*〉,*

_{z}*k,*and

*ɛ*) couple the DPBE equations (Eq. (8)) to the CFD equations (Eqs. (1)–(7)). However, as currently formulated, the CDF equations are solved independent of the DPBE.

*D*

*represents floc diameter of size class*

_{i}*i*(m),

*D*

*is monomer diameter (m),*

_{0}*D*

*is fractal dimension,*

_{f}*k*is lacunarity (generally set as 1, which implies no lacunarity),

*ρ*

*is particle density (kg/m*

_{s}^{3}),

*ρ*

*is fluid density (kg/m*

_{w}^{3}),

*g*is gravitational acceleration (9.81 m/s

^{2}), and

*η*is fluid viscosity (kg/m/s). In Eq. (10),

*2*

*represents the number of monomers forming an*

^{i-1}*i*th particle/floc by following the discretized size classification strategy of the DPBE, which will be described in the following section.

### 2.1.3 Kinetics of particle/floc aggregation and breakage

**(**

*agg/break*

**)**

*). These terms are written as a series of differential equations in Eq. (11). The particle/floc number concentration in a certain discrete size range (*

_{i}*n*

*) is used as a dependent variable of a partial differential equation. Following the discretization scheme of the DPBE, each mean particle size class contains two times the number of constituent monomers in the previous smaller class. Thus if δ is the beginning (irreducible) particle size, class 1 will contain particles of size δ, class 2 will contain particles of size 2δ, class 3 will contain sizes 3δ and 4δ, class 4 will contain 5δ through 8δ, class 5 will contain 9δ through 16δ, and so on. Because the maximum particle size in class*

_{i}*i*increases as 2

^{(i-2)}, 30 classes will contain particles sizes varying from δ to 2

^{28}δ, which represents a growth factor of more than 268 million. Ignoring transport and settling for notational convenience, the partial differential equations with discrete particle/floc size sets may be written as

##### (11)

$$\begin{array}{l}\frac{\partial {n}_{i}}{\partial t}={\left(agg/break\right)}_{i}=\underset{(\text{I})}{\underbrace{{n}_{i-1}\sum _{j=1}^{i-2}{2}^{j-i+1}\alpha (i-1,j)\beta (i-1,j){n}_{j}}}\\ +\underset{(\text{II})}{\underbrace{\frac{1}{2}\alpha (i-1,i-1)\beta (i-1,i-1){n}_{i-1}^{2}}}\\ -\underset{(\text{III})}{\underbrace{{n}_{i}\sum _{j=1}^{i-1}{2}^{j-i}\alpha (i,j)\beta (i,j){n}_{j}}}\\ -\underset{(\text{IV})}{\underbrace{{n}_{i}\sum _{j=i}^{(\text{max}\hspace{0.17em}i)}\alpha (i,j)\beta (i,j){n}_{j}}}-\underset{(\text{V})}{\underbrace{a(i){n}_{i}}}\\ +\underset{(\text{VI})}{\underbrace{\sum _{j=i+1}^{(max\hspace{0.17em}i)+2}b(i,j)a(j){n}_{j}}}\end{array}$$*i*-sized particle/floc generation by collision with other smaller particle/floc classes, (II) generation by collision within the

*i-1*class, (III) disappearance by collision with smaller classes, (IV) disappearance by collision with equal or larger classes, (V) disappearance by fragmentation of the

*i*class, and (VI) generation by fragmentation of larger classes. These mechanisms are illustrated in Fig. 1.

*α, β, a,*and

*b*) are incorporated into the aggregation and break-up kinetics. The collision efficiency factor (

*α*) represents the physicochemical properties of solid and liquid to cause inter-particle attachments (aggregation), whereas the collision frequency factor (

*β*) represents the mechanical properties of fluids to induce inter-particle collisions. In experimental and modeling applications, the collision efficiency factor (

*α*) is generally used as an application-specific fitting parameter, and the collision frequency factor (

*β*) is generally applied as a fixed theoretical function correlated with shear rate (

*β*(

*i, j*) = (

*G/6*)·(

*D*

_{i}*+ D*

*)*

_{j}^{3})[9]. Based on experiments, Ding et al. [36] recently modified these collision efficiency and frequency factors by incorporating the concept of critical size (

*D*

*), which subdivides two different aggregation kinetic regions (the fast and slow aggregation regions) with respect to particle/floc size. In our preliminary modeling and simulation study, application of the critical size concept was found to prevent particle/floc overgrowth beyond the highest particle/floc size and consequently to minimize unexpected mass losses caused by mass escaping out of the defined size range. Thus, in this research, Ding’s critical size (*

_{c}*D*

*) was used as a limiter to prevent unrealistic particle/floc overgrowth in aggregation and break-up kinetics. Eqs. (12) and (13) represent collision efficiency and frequency functions, respectively.*

_{c}*D*

*is the diameter of the*

_{i}*i*th class particle/floc, and

*D*

*is the critical diameter at which 50% of the collisions are successful in forming aggregates [36].*

_{c}##### (12)

$$\alpha (i,j)=\frac{1}{1+{\left(\left({D}_{i}+{D}_{j}\right)/\left(2{D}_{c}\right)\right)}^{3}}$$##### (13)

$$\begin{array}{ll}\beta (i,j)=\frac{1}{6}{\left(\frac{\varepsilon}{v}\right)}^{1/2}{\left({D}_{i}+{D}_{j}\right)}^{3}\hfill & if\hspace{0.17em}{D}_{i},{D}_{j}\le {D}_{c}\hfill \\ \beta (i,j)=\frac{1}{6}{\left(\frac{\varepsilon}{v}\right)}^{1/2}8{\left({D}_{c}\right)}^{3}\hfill & if\hspace{0.17em}{D}_{i},{D}_{j}\ge {D}_{c}\hfill \end{array}$$*a*

*is the selection rate constant, and*

_{0}*V*

*is volume of the*

_{i}*i*th class particle/floc size.

### 2.2. Numerical Simulation

*k–ɛ*turbulence models were selected to simulate flow velocities and turbulence. This resulted in nodal values for (〈

*U*

*〉, 〈*

_{x}*U*

*〉, 〈*

_{y}*U*

*〉,*

_{z}*k,*and

*ɛ*) (Eqs. (1)–(7)). Three different flow conditions that represent low, moderate, and high turbulence conditions were simulated with FLOW-3D, thus data resulting from three steady state flow fields were obtained and saved for the subsequent multi-dimensional DPBE simulation.

*n*

*) of a computational cell was updated at each time step with the inflow and outflow, which are determined by velocities through cell interfaces and concentrations of neighbor computational cells at each time step. Table 1 outlines the numerical scheme used to solve the multi-dimensional DPBE with operator splitting and flux-corrected upwind algorithms. At each time step, particle/floc advection equations were solved with LeVeque flux-corrected upwind scheme and stepwise particle/floc dispersion-reaction equations were calculated implicitly with the Gauss–Seidel iterative method.*

_{i}*Flux*

_{in}*= Flux*

*), while the water surface was treated as a closed boundary (*

_{out}*Flux*

_{out}*= 0*). The bottom layer of the mixing zone was set as a closed boundary for fluid but an open boundary for settling particles/flocs. In other words, for simplification of the model pond system, settling particles/flocs were allowed to move through the bottom layer of the zone, thereby leaving the domain, while fluid remained in the computational domain. The volumetric influent flow rate was initially set to the fixed value of 8 m

^{3}/m/min, which is equivalent to 2.5 min of mean hydraulic residence time (

*t*

_{mean}*= Volume/FlowRate*)within the model mixing zone. However, to create different levels of fluid turbulence and to compare the effects of turbulent intensity on flocculation efficiency, influent flow velocities were set at three different values (0.222, 0.334, and 0.667 m/s) by adjusting inlet width. Influent clay particles (monomers) were assumed to be spheres with a diameter of 1 μm and density of 2.65 g/L. The influent solid concentration was set to 2 g/L, which is equivalent to a particle number concentration of 1.47 × 10

^{15}/m

^{3}.

*D*

*,*

_{f}*D*

*, and*

_{c}*a*

*) were used for aggregation and break-up kinetics. The fractal dimension (*

_{0}*D*

*) was selected as 2.5, which is an intermediate value in the dataset obtained from previous studies [33, 35, 46, 47]. A critical diameter (*

_{f}*D*

*) and breakage kinetic constant (*

_{c}*a*

*) were chosen rather arbitrarily as 100 μm and 10/s following Ding’s recent flocculation theory [36]. However, these constants are previous site-specific values, so it is ultimately recommended that more applicable constants be measured with settling and kinetic experiments appropriate for retention pond applications.*

_{0}### 3. Results and Discussion

*U*〉) and shear rate distributions (

*G*= (

*ɛ*/

*ν*)

^{1/2}), respectively. In the low turbulence condition (case 1), velocity vectors were uniformly directed from the inlet to the outlet, and shear rates were limited to a low level with a maximum shear rate of 13.5/s (Fig. 4(a)). However, in the high turbulence condition (case 3), a swirling zone was identified above the inlet, and high shear rates, with a maximum of 79.3/s, were observed near the inlet (Fig. 4(c)). The moderate turbulent flow condition (case 2) showed flow characteristics intermediate between the two extreme cases (Fig. 4(b)). Later in this paper, we will illustrate the theoretical effects of the velocity and shear rate distributions on flocculation efficiencies.

*Mass*

*,*

_{in,acc}*Mass*

*,*

_{out,acc}*Mass*

*, and*

_{deposit,acc}*Mass*

*represent time-integrated masses caused by in-flux at the inlet, outflux at the outlet, deposition on the bottom, and retention in the pond, respectively, with time progression. Theoretically,*

_{retained}*Mass*

*should be equal to the sum of*

_{in,acc}*Mass*

*,*

_{out,acc}*Mass*

*, and*

_{deposit,acc}*Mass*

*, that is, the mass balance calculated from Eq. (16) should be 100%. Contrary to our expectation, mass balances for low, moderate, and high turbulence conditions were all below, or slightly below, 100% at steady state conditions (99.7%, 97.8%, and 96.1%, respectively). However, these balances were in an acceptable error range, considering the approximating nature and complexity of the numerical methods. In Fig. 5(a), the mass fractions by particle/floc deposition on the bottom (*

_{retained}*Mass*

_{deposit,acc}*/Mass*

*) are also shown for the three different turbulence conditions. The mass fraction by deposition in the high turbulence condition (case 3) was found to be much higher than that in the low turbulence condition (case 1) because high turbulence enhanced flocculation and the subsequent sedimentation processes. These mass fractions and balances became stabilized as the mixing zone systems approached steady state conditions.*

_{in,acc}##### (16)

$$Mass\hspace{0.17em}Balance(\%)=\frac{Mas{s}_{out,acc}+Mas{s}_{deposit,acc}+Mas{s}_{retained}}{Mas{s}_{in,acc}}$$*D*

*), defined by Eq. (17) [48], flowing through the outlet was tracked with time progression to check numerical stability. In Eq. (17)*

_{mm}*m*

*represents the mass of all particles in the*

_{i}*i*th particle/floc size class, and

*M*represents the total mass for all particle/floc size classes.

##### (17)

$${D}_{mm}=\frac{\sum {m}_{i}{D}_{i}}{M}=\left(\frac{{m}_{1}}{M}{D}_{1}+\frac{{m}_{2}}{M}{D}_{2}+\cdots +\frac{{m}_{i}}{M}{D}_{i}\right)$$*D*

*) oscillated and then appeared to stabilize gradually. As mentioned in the previous section, the CFD-DPBE model consists of highly coupled nonlinear equations, which may produce fluctuating results. Strogatz [49] discussed the tendency of such nonlinear equations to produce oscillatory behavior in numerical simulations. A variety of phenomena can contribute to this behavior, including ‘chaotic’ behavior. Thus, the observed oscillatory behaviors shown in Fig. 5(b) were ascribed to the nonlinear nature of the CFD-DPBE model. Such behavior should be examined closely in future experimental and modeling studies.*

_{mm}*D*

*), solid concentration distributions at steady state conditions were investigated in the model-mixing zone. Figs. 6 and 7 show the distributions of mass-weighted mean particle/floc size and solid concentration, respectively, in the three different turbulent flow fields. In case 1, with low turbulence, mass mean particle/floc sizes were less than 27 μm, and solid concentrations were near-homogeneously distributed without particle/floc deposition. However, in case 3 with high turbulence, mass-weighted mean particle/floc sizes grew to 195 μm, which is a sufficient size for the particles/flocs to escape the computational system by settling and being deposited on the bottom. Thus, a longitudinal gradient of solid concentrations was observed in the computational domain due to particle/floc sedimentation. The moderate turbulent flow condition produced results approximately midway between the two extremes. Further, the swirling zones above the inlet in cases 2 and 3 were found to work as small flocculation compartments. Particles/flocs traveling through these swirling zones are more exposed to flocculation and thus tend to grow larger than those passing through the other zones. For example, in case 3, particles/flocs in the swirling zone grew to about 200 μm, while those right next to the swirling zone remained below 50μm.*

_{mm}*D*

*) and deposited mass fraction (*

_{mm}*Mass*

*/*

_{deposit,acc}*Mass*

*) in case 3, with the highest influent flow velocity and shear rate, were up to 7.5 and 12.1 times higher than those in case 1, with the lowest influent flow velocity and shear rate, respectively. As expected, turbulence in sediment retention ponds enhance flocculation efficiency in the mixing zone, at least to a certain extent. In Fig. 8, the cumulative mass distributions of particles/flocs flowing through the outlet are shown for the three different turbulence conditions studied. As expected, for the low turbulence condition (case 1), the particle/floc size distribution was more weighted toward the small size range than that for the moderate and high turbulence conditions (cases 2 and 3). Thus, in case 1, raw clay particles coming through the inlet are not aggregated properly in the turbulent mixing zone, so a large fraction of particles/flocs may not settle appropriately in the subsequent sedimentation basin. In conclusion, considering the results in Table 2 and Fig. 8 from the steady state CFD-DPBE simulations, the turbulence conditions in a turbulent mixing zone were found to have important effects on both flocculation and subsequent sedimentation efficiencies. Optimization of this situation is an important topic for both experimental and theoretical future studies.*

_{in,acc}### 4. Conclusions

The employed CFD software (FLOW-3D) was a useful tool for generating steady state flow field data, such as flow velocities and shear rates, which were used in subsequent multi-dimensional DPBE simulations.

As an alternative to QMOM, DPBE formulation was applied to simulate a multi-dimensional flocculation/sedimentation process. Solution of the multi-dimensional DPBE provides more readily understandable results for engineers and scientists with slightly more computations than QMOM but well within the capabilities of modern personal computers for two-dimensional flow fields.

A standard, central-differencing, finite difference approach was judged inadequate for simulating the flocculation and sedimentation processes in sediment retention ponds due to computational instability caused by nonlinearity, advection dominance, and complexity of the DPBE model. Thus, operator-splitting and LeVeque flux-corrected algorithms were applied to overcome the computational difficulties. The detailed numerical model is available from the authors upon request.

In application of the CFD-DPBE model, increased turbulence was found to enhance flocculation and sedimentation efficiencies. However, methodology for optimizing this effect requires further study.