**Experimental Investigation of Erosion Characteristics of Fine-Grained Cohesive Sediments**

**Bommanna Gounder Krishnappan 1,\*, Mike Stone 2, Steven J. Granger 3, Hari Ram Upadhayay 3, Qiang Tang 3, Yusheng Zhang <sup>3</sup> and Adrian L. Collins <sup>3</sup>**


Received: 28 April 2020; Accepted: 22 May 2020; Published: 25 May 2020

**Abstract:** In this short communication, the erosion process of the fine, cohesive sediment collected from the upper River Taw in South West England was studied in a rotating annular flume located in the National Water Research Institute in Burlington, Ontario, Canada. This study is part of a research project that is underway to model the transport of fine sediment and the associated nutrients in that river system. The erosion experimental data show that the critical shear stress for erosion of the upper River Taw sediment is about 0.09 Pa and it did not depend on the age of sediment deposit. The eroded sediment was transported in a flocculated form and the agent of flocculation for the upper River Taw sediment may be due to the presence of fibrils from microorganisms and organic material in the system. The experimental data were analysed using a curve fitting approach of Krone and a mathematical model of cohesive sediment transport in rotating circular flumes developed by Krishnappan. The modelled and measured data were in good agreement. An evaluation of the physical significance of Krone's fitting coefficients is presented. Variability of the fitting coefficients as a function of bed shear stress and age of sediment deposit indicate the key role these two factors play in the erosion process of fluvial cohesive sediment.

**Keywords:** erosion; cohesive sediments; rotating circular flume; mathematical modelling; fitting coefficients; sediment deposition; flocculation; bed shear stress; consolidation

#### **1. Introduction**

Fine-grained cohesive sediment plays an important role in the transportation of pollutants and it is a key driver of water quality degradation in rivers. Rigorous quantification of cohesive sediment transport processes is fundamental for predicting sediment and associated contaminant transport in aquatic systems (Horowitz and Elrick [1]; Luoma and Rainbow [2]). The development of reliable numerical models to simulate cohesive sediment transport dynamics requires an accurate description of fundamental sediment transport processes such as erosion, deposition, and transport of solids in suspension (Grabowski et al. [3]). Factors affecting the erosion characteristics of cohesive sediments include the rate of bed shear (Amos et al. [4]), the degree of consolidation/age of deposit (Lick and McNeil [5]), bio-stabilisation effects by microorganisms (Friend et al. [6]) and the initial conditions that created the deposit (Lau et al. [7]). In a re-examination of the sediment erosion data from flume studies conducted by Roberts et al. [8] and Zreik et al. [9], Krone [10] demonstrated the importance of the sediment bed matrix structure on erosion characteristics of fine sediment deposits. Factors controlling the deposition characteristics of cohesive sediment include the phyco-chemical properties of sediment water mixture, the concentration of organic matter including microorganisms in the water column and the turbulence characteristics of the flow field which in turn, can cause the suspended sediment particles to flocculate and change their porosity, density and settling characteristics. At the present state of knowledge, numerical models of cohesive sediment transport rely mainly on laboratory experiments using specialised flumes such as a rotating annular flume for the determination of transport parameters that include the critical shear stresses for erosion and deposition, erosion and deposition rates and properties of sediment flocs for site-specific sediments.

A research project is underway in the upper River Taw in South West England to model the transport of fine-grained cohesive sediments. The upper River Taw has been instrumented as a landscape scale observatory for exploring the interactions between climate, land use, farm management and water quality in conjunction with a larger strategic programme exploring pathways for improving the sustainability of agriculture. The upper River Taw drains both organic-rich peaty and podzolic upland moorland soils near its source and clayey gley and brown earth soils in the more intensive agricultural areas (beef and sheep grazing, cereals, maize, oil seed rape) on the more intensively farmed lowland adjacent to the moor. Sediment stress in the study area has been highlighted as a critical factor impacting on lithophilous fish species dependent upon clean bed gravels for their early life stages including the incubation of progeny in redds cut into bed gravels.

As part of the research project on modelling cohesive sediment transport, a survey of 18 cross-sections located ~0.5 km apart was conducted in the upper River Taw (Figure 1) during the summer low flows in 2018. The bank-full width of the study river is ~10 m and the average depth is ~1 m. The riverbed is armoured with large stones including pebbles and cobbles. Bedrock outcrops are also observed. Although the matrix bed material in the study reach is very coarse, fine-grained sediment is present in the river channel bed because of the entrapment process, which retains fine material in the lee of large rocks or directly within the gravel bed matrix as interstitial fines. Fine-grained sediment deposits are also present in recirculation eddy zones along the edges of the river channel where the bed shear stresses of the flow field are low. Representative samples of fine-grained sediment from the study reach were collected using a conventional bed sediment remobilisation technique (Lambert and Walling [11]; Duerdoth et al. [12]; Naden et al. [13]) at all 18 cross sections and 18 one L bottles of sediment slurry were then shipped to Canada for testing in a rotating annular flume located in the National Water Research Institute in Burlington, Ontario, Canada.

In this short communication, preliminary results from this research project on cohesive sediment transport are presented. Specifically, a component of the cohesive sediment transport, namely, the erosion process, was investigated by carrying out erosion experiments using a rotating annular flume. A novel approach was used to analyse the erosion data. The approach is based on a methodology proposed by Krone [10] and a mathematical model of cohesive sediment transport in rotating annular flumes developed by Krishnappan [14].

**Figure 1.** Map of the upper River Taw study catchment, showing location in south west England, channel bed cohesive sediment sampling locations and flow gauging stations.

#### **2. Material and Methods**

#### *2.1. Rotating Annular Flume*

The rotating circular flume used in this study was 5.0 m in mean diameter, 30 cm wide and 30 cm deep and it rests on a rotating platform which was 7 m in diameter. An annular lid fitted inside the flume with close tolerance (about 1 mm gap all around) and it rotated in the opposite direction to the flume's rotation. The annular lid maintained contact with the water surface within the flume during experiments. The generated flow fields in such assemblies were nearly two dimensional with near constant bed shear stress across the channel and with minimum secondary circulation in the transverse direction (Petersen and Krishnappan [15]). The flow field in this rotating annular flume assembly was computed using the 3D hydrodynamic flow model, PHOENICS (Rolston and Spalding [16]). The bed shear stress computed by the model was verified using direct measurements of bed shear stress using a preston tube (Krishnappan and Engel [17]). The main advantage of rotating flumes over linear flumes is that detrimental effects of the pump and the pipe system on the floc structure of the cohesive sediment were avoided, thereby permitting reliable studies on floc behaviour to be conducted. A complete description of the flume can be found in Krishnappan [18].

The sediment samples shipped from the UK were stored in a cold room until the start of the experiments. A large composite sediment sample was prepared by combining all 18 bottles of samples and wet sieved using a 200 mesh sieve before being added to the flume. The erosion experiments were conducted by using the standard procedure which involves operating the flume at high speed initially to suspend the sediment and then lowering the flume speed gradually to allow the suspended sediment to deposit and form a uniform bed. The bed was then allowed to age for three different time periods (herein referred to as age of deposit), namely, 22 h, 38 h and 160 h. During an erosion test, bed shear stresses were increased over time in steps as a stair-case function. The shear stress steps used were 0.06, 0.09, 0.12, 0.17, 0.27 and 0.33 Pa, with each maintained for a period of one hour. In each shear stress step, sediment samples were collected every 10 min to determine the concentration of the eroded sediment as a function of time. When the concentration reached a steady-state value, the flume speed was increased to the next level. Whenever there was sufficient sediment suspended in the water column, sediment samples were collected for size analysis using a LISST (Laser In-situ Scattering and Transmissometry 100X; Sequoia Scientific, Bellevue, WA, USA) particle size analyser and an image analysis system. This procedure was repeated until the maximum permissible flume speed was reached.

#### *2.2. Methodology of Krone*

Krone [10] developed his methodology when he analysed erosion data reported by Zreik et al. [9] who studied the erosion behaviour of cohesive sediments from Boston Harbour using a rotating circular flume. Measured sediment concentration data from Zreik et al. [9] showed that the concentration of eroded sediment increased rapidly at the beginning of each applied shear stress step change but then levelled off with time before the next incremental step in shear stress.

Krone [10] fitted the following equation to describe the variation in sediment concentration during a shear stress step:

$$C = \frac{\frac{1}{c\_1}t}{\frac{c\_0}{c\_1}t + 1} \tag{1}$$

where *C* = incremental concentration of eroded sediment during a shear stress step; t = time after the shear stress change; *c*<sup>0</sup> and *c*<sup>1</sup> = constants. This equation fitted the experimental data of Zreik et al. [9] very well and, importantly, the constants *c*<sup>0</sup> and *c*<sup>1</sup> have physical meaning. For example, when t tends to infinity, the ratio 1/*c*<sup>0</sup> takes on the value of the concentration near the end of the shear stress step in question. When t = 0, the ratio 1/*c*<sup>1</sup> assumes the value of the erosion rate, i.e., *dC*/*dt* = 1/*c*<sup>1</sup> at *t* = 0. The constants *c*<sup>0</sup> and *c*<sup>1</sup> were evaluated by knowing sediment concentrations at two different times (near the beginning and near the end of the shear stress step) using the following relationships:

$$c\_1 = \left(\frac{1}{C\_T} - \frac{1}{C\_L}\right)T \qquad\qquad\text{and } c\_0 = \frac{1}{C\_L} \tag{2}$$

where *CT* = concentration near the beginning of the shear stress step, *CL* = limit concentration at the end of the shear stress step and *T* is the elapsed time from the beginning of the shear stress step till the time when *CT* is specified. Knowing *c*<sup>0</sup> and *c*1, Equation (1) can be used to calculate the concentration of the eroded sediment for the entire duration of each shear stress step. In addition, Equation (1) is used to derive an erosional rate function that was used in the model of cohesive sediment transport developed by Krishnappan [2]. The erosion rate of sediment can be expressed as follows:

$$E = h \frac{d\mathbf{C}}{dt} \tag{3}$$

where *E* is the erosion rate in gm/m2s, h is the depth of water in the flume in metres and *dC*/*dt* is the concentration gradient which can be evaluated from Equation (1). Substituting the expression of *dC*/*dt* in Equation (3), the erosion rate function can be derived as:

$$E = h \frac{\left(1/c\_1\right)}{\left(\frac{c\_0}{c\_1}t + 1\right)^2} \tag{4}$$

#### *2.3. Mathematical Model of Cohesive Sediment Transport Developed by Krishnappan*

A mathematical model of cohesive sediment transport in the rotating circular flume (called the FLUME model) was used to simulate the erosion process of the sampled upper Taw River sediment. The FLUME model incorporates the erosion rate function of Krone [10] i.e., Equation (4). A full description of the model can be found in Krishnappan [14]. Here, a summary of the salient features of the model is outlined for the sake of completeness. The FLUME model treats the motion of sediment in the rotating flume in two stages: a transport/settling stage and a flocculation stage. Some salient features of these two fundamental stages of sediment transport are briefly discussed below:

#### 2.3.1. Settling Stage

The governing equation for the settling stage was obtained from mass balance considerations. Assuming that the flow in the rotating flume is uniform in the longitudinal and tangential directions, the one-dimensional mass balance equation as shown below was adopted:

$$\frac{\partial \mathbf{C}\_k}{\partial t} + w\_k \frac{\partial \mathbf{C}\_k}{\partial z} = \frac{\partial}{\partial z} \Big(\Gamma \frac{\partial \mathbf{C}\_k}{\partial z}\Big) + S \tag{5}$$

where, *Ck* is the concentration of sediment in size fraction *k*, *wk* is the settling velocity of that fraction and Γ is the dispersion coefficient in the vertical direction. The symbol *t* represents the time axis and *z* represents the coordinate axis in the vertical direction. *S* is the source/sink term. The boundary conditions used for the settling stage are:

At the free surface:

$$-w\_k \mathbb{C}\_k - \Gamma \frac{\partial \mathbb{C}\_k}{\partial z} = 0 \tag{6}$$

At the bed:

$$-w\_k \mathbb{C}\_k - \Gamma \frac{\partial \mathbb{C}\_k}{\partial z} = q\_\ell + q\_d + q\_{entrap} \tag{7}$$

The free surface boundary condition expresses the balance between the settling flux and the dispersive flux so that there is no external input of sediment at this boundary. At the bed surface, it is assumed that the settling and dispersive fluxes are balanced by the net amount of sediment exchanged at the sediment-water interface. The sediment exchange at the sediment-water interface can occur when: (1) sediment is eroded from the bed and entrained into the water column (*qe*) (2) sediment settles to the bed and stays on the bed as the deposited sediment (*qd*) (3) sediment settles to the bed and a portion of the deposited sediment can ingress into the interstitial pores of the gravel bed and become unavailable for further erosion and entrainment (*qentrap*).

The first two components of the sediment exchange at the sediment-water interface have been studied extensively in the literature, including by Partheniades [19], Mehta and Partheniades [20], Parchure [21], Lick [22], Krone [10], Krishnappan [23], Krone [24], Mehta and Partheniades [20], and Lick [22]. In this study, the approaches proposed by Krone [10] and Krone [24] for erosion and deposition respectively have been adopted. Accordingly, the erosion flux was calculated as:

$$q\_c = E \tag{8}$$

where the erosion rate, *E* is given by Equation (4) and the deposition flux was calculated as:

$$q\_d = P w\_k C\_{kb} \tag{9}$$

where *P* is a probability parameter which gives a measure of the probability that a sediment particle settling to the bed, stays at the bed. *Ckb* is the near-bed concentration of the sediment fraction *k*. Krone [24] proposed a relationship for *P* as:

$$P = \left(1 - \frac{\pi}{\tau\_{crd}}\right) \tag{10}$$

where τ is the bed shear stress and τ*crd* is the critical shear stress for deposition, which is defined as the bed shear stress above which none of the initially suspended sediment would deposit.

The third component, namely, the entrapment component, is assumed to be directly proportional to the settling flux near the bed and is expressed as:

$$
\mathfrak{q}\_{entmp} = \alpha \mathfrak{z} \mathfrak{w}\_{k} \mathbb{C}\_{k\flat} \tag{11}
$$

where α is the proportionality constant and is termed the entrapment coefficient. The entrapment coefficient is expected to be a function of porosity of the gravel substrate, the thickness of the gravel bed layer and the permeability of the substrate. Further research is needed to quantify this term as a function of the bulk properties governing the entrapment process.

#### 2.3.2. Flocculation Stage

The flocculation stage was modelled using the coagulation equation (Fuchs [25]) shown below:

$$\frac{\partial \mathcal{N}(i,j)}{\partial t} = -\beta \mathcal{N}(i,t) \sum \prescript{\circ}{}{\imath}\_{i-1} \mathcal{K}(i,j) \mathcal{N}(j,t) + \frac{1}{2} \beta \sum \prescript{\circ}{}{\imath}\_{i-1} \mathcal{K}(i-j,j) \mathcal{N}(i-j,j) \mathcal{N}(j,t) \tag{12}$$

where *N*(*j*, *t*) and *N*(*j*, *t*) are number concentrations of particle size classes *i* and *j*, respectively. *K*(*i*, *j*) is the collision frequency function, which is a measure of the probability that a particle of size *i* collides with a particle of size *j* in unit time and β is the cohesion factor, which defines the probability that a pair of collided particles will coalesce and form a new floc. The β term accounts for the effects of cohesion properties such as the electrochemical properties of the sediment-water mixture and the effects of polymers secreted by microorganisms. In this study, is treated as an empirical parameter and was determined through the calibration of the model using experimental data from the flume.

The first term on the right side of Equation (12) gives the reduction of particles in size class *i* because of flocculation of these particles in size class *i* with all other size classes. The second term gives the generation of new particles in the size class *i* because of flocculation of particles in smaller size classes. In this process, the mass of particles is assumed to be conserved. Equation (12) was simplified by considering the particle sizes in discrete size classes where the continuous particle size space is considered in discrete size ranges. Each range is considered as a bin containing particles of a certain size range. For example, *ri* is the geometric mean radius of particles in bin 1. The particle ranges were selected in such a way that the mean volume of particles in bin *i* is twice that of the preceding bin.

Under this scheme, when particles of bin *i* flocculate with particles in bin *j* (*j* < *i*), the newly formed particles will fit into bins *i* and *i* + 1. The proportion of particles going into these bins can be calculated by considering the mass balance of particles before and after flocculation. With this simplification, the coagulation equation can be expressed in discrete form as follows:

$$\frac{\Delta N\_i}{\Delta t} = -\sum \beta \mathbf{K}(r\_i, r\_\rangle) \mathbf{N}\_i \mathbf{N}\_\rangle + \sum f\_{i,\dagger} \beta \mathbf{K} \{ r\_{i\prime} r\_\rangle \} \mathbf{N}\_i \mathbf{N}\_\rangle + \sum (1 - f\_{i,\dagger}) \beta \mathbf{K} \{ r\_{i\prime} r\_\rangle \} \mathbf{N}\_{i-1} \mathbf{N}\_\rangle \tag{13}$$

where *fi*,*<sup>j</sup>* is the allocation function given by:

$$f\_{i,j} = \left(\rho\_{s,j}V\_i + \rho\_{s,j}V\_j - \rho\_{s,j+1}V\_{i+1}\right) / \left(\rho\_{s,i}V\_i - \rho\_{s,j+1}V\_{i+1}\right) \tag{14}$$

where, ρ*<sup>s</sup>* denotes the density of the flocs and *V* is the volume of the flocs. The density of the flocs is dependent on floc size and, here, an empirical relationship proposed by Lau and Krishnappan [26] was adopted in this study. The form of the relationship is as follows:

$$
\rho\_{\mathbb{S}} - \rho = \left(\rho\_{\mathbb{P}} - \rho\right) \exp(-bD^{\mathbb{C}}) \tag{15}
$$

where ρ*s*, ρ and ρ*<sup>p</sup>* are densities of flocs, water and parent material forming the flocs, respectively. *D* is the diameter of the floc and b and c are empirical coefficients to be determined through model calibration. The settling velocity of flocs, needed for the settling stage of the sediment particle motion, was calculated using the above density relation and Stoke's Law. The resulting expression is as follows:

$$w\_k = \left(\frac{1.65}{18}\right) \left(\frac{gD\_k^2}{\upsilon}\right) \exp\left(-bD^c\right) \tag{16}$$

where *g* stands for acceleration due to gravity and υ stands for the kinematic viscosity of water.

The collision frequency function *K*(*i*, *j*) takes different functional forms depending on the collision mechanisms that bring the particles to close proximity. The mechanisms considered in this study include (1) Brownian motion (*Kb*); (2) turbulent fluid shear (*Ke*h); (3) inertia of particles in turbulent flow (*KI*), and; (4) differential settling of flocs (*Kds*). An effective collision frequency function, *Ke f* , was calculated in terms of the individual collision functions as follows:

$$K\_{\varepsilon f} = K\_b + \sqrt{\left(K\_{sh}^2 + K\_l^2 + K\_{ds}^2\right)} \tag{17}$$

The collision frequency functions for the individual collision mechanisms can be found in Valioulis and List [27].

The break-up of flocs due to turbulent fluctuations of the flow was modelled using a scheme proposed by Tambo and Watanabe [28]. According to this scheme, a "collision-agglomeration" function was introduced as a multiplier to the collision frequency function *Ke f* . This produced an effective collision frequency that resulted in an optimum floc size distribution for a given level of turbulence. The function proposed by Tambo and Watanabe [28] is as follows:

$$
\alpha\_{\mathcal{R}} = \alpha\_0 \left( 1 - \frac{R}{R\_{\mathcal{M}} + 1} \right)^n \tag{18}
$$

where *R* is the number of primary particles contained in a floc and *Rm* is the number of particles contained in the biggest floc for the given level of turbulence. The parameters α<sup>0</sup> and *n* assume values of 1/3 and 6, respectively. The flow field and the turbulence characteristics in the rotating circular flume needed for the FLUME model were predicted using PHOENICS (Rosten and Spalding [16]).

The FLUME model was used in this study to simulate the erosion process, even though the model is capable of predicting the total transport of cohesive sediment in rotating circular flumes including the deposition and the entrapment processes. The FLUME model, therefore, can be used to determine all parameters governing the transport of cohesive sediments by carrying out experiments using a rotating circular flume and then use these parameters in models that can predict the transport of cohesive sediment under field conditions.

#### **3. Results and Discussion**

The results of erosion tests from all three runs carried out for the present study are presented in Figure 2. The erosion data indicate that the sediment deposit was completely stable until the shear stress reached a value of 0.09 Pa, which can be considered as the critical shear stress for the erosion of the surficial layer of the sediment deposit. The data presented also show that six out of the seven shear stress steps at each age of deposit produced erosion of the riverbed sediment. At each of the six shear stress steps when erosion of the sediment bed occurred, the eroded sediment concentration increased suddenly after the shear stress change, suggesting that the sediment erosion follows the pattern of "bulk" erosion. As the sediment erosion continues, the erosion rate decreases due to the increasing strength of the bed with bed depth. Erosion during the later stages of the shear stress step is defined as the "surface" erosion when the removal of the material from the bed is particle by particle and the concentration in the water column approaches a steady state value. Additionally, the influence of the age of deposit is evident for the lower shear stress steps, since the concentration of the eroded sediment decreased as the age of deposit increased (Figure 2). For higher shear stress steps (i.e., for 0.27 and 0.33 Pa), the influence of age of deposit is not as pronounced as for the lower shear stress steps.

**Figure 2.** Results from the erosion tests.

The size distribution of the eroded sediment measured using a LISST and the image analysis system indicate that the sediment is flocculated. Figure 3 shows the size distribution data measured using LISST for the shear stress step of 0.33 Pa and where the age of deposit is 160 h.

**Figure 3.** Size distribution of eroded sediment at the shear stress step of 0.33 Pa and age of deposit equal to 160 h.

The undisturbed and sonicated size distribution of eroded sediment for one set of experimental conditions (0.33Pa and deposit age of 160 h) are presented in Figure 3. The dotted line represents the size distribution of the eroded sediment as determined by the LISST for the undisturbed sample, whereas the solid line represents the distribution for the sonicated sample. Sonication breaks up all the flocs that are present in the sample. From these two distributions, we can conclude that the eroded sediment in the water column is in a flocculated state and the sediment exhibits cohesive tendencies. The change in the predominant size of the flocculated sediment from ~50 μm to ~25 μm after sonication highlights how flocculation can influence particle morphology and, by extension, volumetric concentration (Figure 3). A particle size of 63 microns is considered in the literature as the division between cohesive and non-cohesive sediments.

A photomicrograph of the eroded sediment sample collected for the image analysis is shown in Figure 4.

**Figure 4.** A photomicrograph of the eroded sediment for shear stress step of 0.33 Pa and age of deposit equal to 160 h.

This figure provides visual evidence of flocculation of the particles and also shows that mico-flocs are the building blocks of larger flocs as previously reported by Stone et al. [29]. The particles are interconnected through loose fibril material that might have been secreted by the microorganisms or the organic material that is present in the system. Therefore, flocculation is enhanced by microorganisms and the organic material and this process will have to be accurately represented in models to simulate sediment transport in the upper River Taw River.

Krone's curve fitting approach was used for the data shown in Figure 2. For each shear stress step where the sediment erosion was present, the constants *c*<sup>0</sup> and *c*<sup>1</sup> were calculated using Equation (2) with T = 10 min. The values of these constants, termed herein as "fitting coefficients" are listed in Table 1.


**Table 1.** Fitting coefficients *c*<sup>0</sup> and *c*<sup>1</sup>

With values of *c*<sup>0</sup> and *c*1, Equation (1) was applied to all the shear stress steps and the concentration variation as a function of time was calculated. A comparison of the fitted concentration variation with the measured data is presented in Figure 5 for a shear stress step of 0.17 Pa and consolidation time of 38 h as an example. The modelled and measured data are in agreement for all six shear stress steps in all three sediment consolidation time tests.

**Figure 5.** Comparison between measured and fitted sediment concentrations using Equation (1) for shear stress step: 0.17 Pa; Age of deposit: 38 h (r2 = 0.993).

The FLUME model was applied to the present erosion experiments to simulate the concentration of the eroded sediment in the water column for all three experiments. The model was applied to each shear stress step using the appropriate parameters, *c*<sup>0</sup> and *c*1, corresponding to that particular shear stress step and the age of the deposit. The deposition and entrapment fluxes were suppressed for the present simulations. The simulated concentration of the eroded sediment is compared to the measured concentration in Figure 6 and a favourable agreement between the two was observed. The erosion rate function that was derived from the fitting equation (Equation (1)) of Krone [10] works well with the model able to produce sediment concentrations that agree well with the measured data.

**Figure 6.** Comparison of measured and modelled sediment concentrations predicted by the FLUME model.

The fitting coefficients, c0 and c1, have different values for each shear stress step and age of deposit (Table 1). The variability of these coefficients with each shear stress and age of deposit is shown in Figure 7 which shows that these coefficients vary as a function of both shear stress and age of the deposit.

The coefficients exhibit a complex behaviour with respect to bed shear stress. Both coefficients decrease initially as the shear stress increases, reach minimum values and then increase as the shear stress is increased further. Notably, this behaviour was observed for all three ages of deposit tested. The minimum conditions for the coefficients imply maximum erosion rate and the maximum amount of sediment eroded (implying "bulk" erosion). The minimum condition for the coefficients shifts to the right as the age of deposit increases (Figure 7). Notably, when the age of deposit is 22 h, minimum conditions for both coefficients are at around 0.15 Pa, whereas for the age of deposit of 160 h, the minimum conditions are shifted to 0.21 Pa for *c*<sup>0</sup> and to 0.25 Pa for *c*1. Accordingly, this suggests that the bed shear strength has increased due to the age of the deposit and higher shear stresses are needed to cause "bulk" erosion. The variability of the fitting coefficients demonstrates the importance of the roles that the bed shear stress and the age of deposit play in determining the erosion behaviour of the upper River Taw sediment. This finding will be useful in developing a modelling framework for predicting cohesive sediment transport in the upper River Taw.

**Figure 7.** The variability of the fitting constants as a function of shear stress and age of deposit (**a**): 22 h; (**b**): 38 h and (**c**): 160 h.

#### **4. Conclusions**

A better understanding of cohesive sediment transport processes is required to develop models that can simulate scenarios that address both scientific and policy questions. In this short communication, the erosion process of the fine, cohesive sediment collected from the upper River Taw in South West England was studied by conducting erosion experiments using a rotating annular flume located in the National Water Research Institute in Burlington, ON, Canada. The erosion experimental data show that the critical shear stress for erosion of the upper River Taw sediment is about 0.09 Pa and it did not depend on the age of the deposit. The size distribution measurements for the eroded sediment indicates that the eroded sediment was transported in a flocculated form and hence the upper River Taw sediment exhibits cohesive properties. Photomicrographs of the eroded sediment samples obtained by image analysis suggest that the agent of flocculation for the sampled sediment may be due to the presence of fibrils from microorganisms and organic material present in the River Taw system. The experimental data from the erosion experiments were analysed using an approach proposed by Krone [10] and a mathematical model of cohesive sediment transport in rotating annular flumes developed by Krishnappan [14]. The approach of Krone [10] was applied to the experimental data and the fitting coefficients were established as a function of bed shear stress and age of sediment deposit. Using the approach of Krone [1], an erosion rate function was calculated and was used in the FLUME model of Krishnappan [14]. A comparison of the model predictions and the experimental data confirmed that the fine sediment transport model is capable of simulating accurately the erosion experiments in a rotating circular flume. The variability of the fitting coefficients, as a function of bed shear stress and the age of the deposit, was also examined in the present study. Future work will incorporate the experimental results reported herein into a catchment scale cohesive sediment transport modelling framework for exploring the implications of future climate and land management scenarios on sediment transport and behaviour.

**Author Contributions:** Conceptualization, B.G.K., M.S. and A.L.C.; methodology, B.G.K.; Software, B.G.K.; Validation, S.J.G., H.R.U., Q.T. and Y.Z.; formal analysis, B.G.K.; investigation, S.J.G., H.R.U., Q.T. and Y.Z.; resources, S.J.G.; data curation, S.J.G., H.R.U., Q.T. and Y.Z.; writing-original draft preparation, B.G.K.; review and editing, M.S. and A.L.C.; visualization, S.J.G.; supervision, M.S. and A.L.C.; project administration, M.S. and A.L.C.; funding acquisition, A.L.C. and M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** Rothamsted Research received strategic funding from UKRI-BBSRC (UK Research and Innovation—Biotechnology and Biological Sciences Research Council) and the contribution of institute staff to this paper was funded by the Soil Nutrition institute strategic programme project 3 grant (BBS/E/C/000I0330). Funding for the flume experiments was from an NSERC Discovery Grant RGPIN-2020-06963 awarded to M. Stone.

**Acknowledgments:** The authors acknowledge the cooperation of Ian Droppo of Environment and Climate Change Canada and the use of the Rotating Circular Flume for their study. The flume experiments were carried out by Robert Stephens, a retired Research Technician of Environment and Climate Change Canada.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Integrated Hydrological-CFD Model for Estimating Bacterial Levels in Stormwater Ponds**

#### **Farzam Allafchi 1, Caterina Valeo 1,\*, Jianxun He <sup>2</sup> and Norman F. Neumann <sup>3</sup>**


Received: 1 May 2019; Accepted: 10 May 2019; Published: 15 May 2019

**Abstract:** A hydrological model was integrated with a computational fluid dynamics (CFD) model to determine bacteria levels distributed throughout the Inverness stormwater pond in Calgary, Alberta. The Soil Conservation Service (SCS) curve number model was used as the basis of the hydrological model to generate flow rates from the watershed draining into the pond. These flow rates were then used as input for the CFD model simulations that solved the Reynolds-Averaged Navier-Stokes (RANS) equations with *k-*ε turbulence model. *E. coli*, the most commonly used fecal indicator bacteria for water quality research, was represented in the model by passive scalars with different decay rates for free bacteria and attached bacteria. Results show good agreement with measured data in each stage of the simulations. The middle of the west wing of the pond was found to be the best spot for extracting water for reuse because it had the lowest level of bacteria both during and after storm events. In addition, only one of the four sediment forebays was found efficient in trapping bacteria.

**Keywords:** stormwater reuse; SCS curve number; CFD; fecal indicator bacteria; *E. coli*

#### **1. Introduction**

In recent decades stormwater has been considered as an alternative water source for reuse, specifically for applications that need less than pristine water quality. Reusing stormwater is more critical in water-scarce regions and regions where rainfall patterns and rainfall frequencies are changing [1]. Several regions are trying to reuse stormwater as a sustainable method of water resources management; thus, prompting research into the feasibility of reusing stormwater from a water quality perspective [2].

Although stormwater ponds are built with the primary objective of reducing runoff quantities in order to protect urban areas against flooding, they also improve the quality of stormwater as well [1]. The stormwater quality within a pond varies both spatially and temporally and is not only a function of the quality of the influent but also a function of local hydrological conditions and on the pond's design [3,4]. A water quality study of a very large stormwater pond in Calgary, Alberta, Canada, showed that the quality of stormwater in the pond observed over a three year period did meet the irrigation water quality requirements but only under certain circumstances [2]. Thus, stormwater recycling with pond water often requires continuous or intermittent water quality monitoring of the pond water in order to remain compliant with local regulations. Highly distributed water quality sampling in stormwater ponds is often impractical due to the sizes of these ponds and the cost. In addition, most ponds are not designed with reuse in mind and the extraction point is often located in an ad hoc fashion and possibly in a region of the pond which has higher pollution levels relative to

the rest of the pond because of local hydrodynamic conditions. If retrofitting a pond with the intent to recycle the water, the municipality could undertake a water quality sampling program that collected samples distributed throughout the pond over a period of time in order to identify the optimum location for extracting the "cleanest" water in the pond (assuming there is no treatment of this water). A more cost-effective alternative is to develop a physical model to estimate the bacteria level in the pond that incorporates the factors leading to bacterial contamination of stormwater in retention ponds.

Water quality is impaired when mobilized sources of contaminants are transported away from their original location by runoff and discharged in aquatic environments [5,6]. Fine particles are either detached from the soil or washed off an impervious surface [7]. Determining levels of pathogenic microorganisms is expensive and difficult due to their large diversity; thus, microbiological water quality monitoring procedures often use fecal indicator bacteria (FIB). FIB are present in feces of human and warm-blooded animals in large numbers and can be easily detected [8]. Total coliforms (TC) and fecal coliforms (FC) were considered the main group of FIB during the twentieth century. However, nowadays, *Escherichia coli* (*E. coli*) and intestinal enterococci (IE) are enumerated and considered as FIB because some water-related epidemiological studies have shown that they are better bacterial indicators for predicting sanitary risk [9,10].

The aim of microbiological water quality monitoring of stormwater ponds is primarily to assess the level of fecal contamination in an aquatic system. This contamination varies spatially and temporally throughout the pond and, therefore, very large data sets must be collected to enable an adequate understanding of contamination levels and the processes leading to the degree of pollution in the system. Collecting such a large data set is often intractable, and thus models are often the only feasible approach for gaining greater insight into aquatic ecosystem processes. Two main types of models have been proposed to estimate bacterial levels in aquatic environments: data-driven models, which are also considered black box models, and process-based models [11].

Data-driven models use statistical methods or computational intelligence and machine learning to relate the involved parameters to the state variables (input, internal and output variables) with only a limited number of assumptions about the physical behavior of the system. In contrast, process-based methods predict the FIB concentration based on the mathematical description of sources, sinks, and internal processes influencing FIB levels. The fundamental principle underlying process-based models is conservational equations. The most important factors influencing microorganism fate in the aquatic environment are (i) environmental parameters, and (ii) whether or not they are attached to particles [11]. Few estimates of the percentage of particles with attached bacteria in any collected stormwater sample can be found in the literature. One study found up to 10–28% and 22–30% of fecal coliform and *E. coli*, respectively, were associated with suspended solids in stormflow samples at the mouth of a canal discharging water to a lake with brackish water [12]. Another study found that 30–55% of both fecal coliform and *E. coli* were attached to sediment particles in stormwater samples [13].

The fate and transport of bacteria depend heavily on the attachment to suspended solids. Attachment to sediment protects bacteria from some processes that may accelerate death and decay, such as sunlight and predation. Thus, any modelling should use different decay rates for bacteria attached to sediment vs. free-floating bacteria. For example, the decay rate of free-floating *E. coli* was considered to be twice that of *E. coli* attached to sediment in a study in the Scheldt drainage network in Belgium [14]. In another study in the Blackstone River watershed in Massachusetts, the ratio of the decay rate of free-floating bacteria to attached bacteria was considered to be 4 [15]. Some studies even assumed that attached pathogens did not decay at all [11].

The fate of bacteria is most commonly modeled by first-order kinetic decay, and this first-order model was used to study the efficiency of waste stabilization ponds. Shilton [16] used a pulse tracer study incorporated with computational fluid dynamics (CFD) to find the retention time of a waste stabilization pond in order to use the first-order kinetic decay based on the retention time [16]. In another study on an anaerobic lagoon, decay was calculated based on the retention time resulting from CFD simulations [17]. In both the study of the waste stabilization pond [16] and the anaerobic

lagoon, the objective was only to study the inlet and outlet bacteria concentrations, and they neglected to consider variation within the body. In addition, in both of the studies, the simulation was assumed to be steady-state. Therefore, the retention time, which did not change with time could be an indicator of bacteria concentration. Two other studies on waste stabilization ponds used tracers integrated with CFD simulations to study different configurations of baffles [18,19]. They implemented first-order kinetic decay by a source term in the transport equation of the tracer. All of the abovementioned studies simulated the flow as steady-state and neglected wind. The fluid flow and bacteria fate and transport in stormwater ponds are intrinsically unsteady. Also, the wind cannot be neglected. In addition, the main goal of the present study is to study bacteria concentration within stormwater ponds for the purposes of determining the optimum point at which water may be withdrawn for reuse—the optimum point is that location with the cleanest water and having the lowest level of bacteria concentrations.

In the present study, a comprehensive model was developed that integrates a hydrological model for a catchment draining to a stormwater pond, with a computational fluid dynamics (CFD) model simulating the pond's hydrodynamics. The results of the hydrological model were used as inputs to the CFD simulation. The model results are validated against data collected at the pond. The overall goal of this paper was to enhance knowledge of bacteria fate and transport in stormwater ponds. Furthermore, the developed modeling approach leading to this enhanced understanding may ultimately be used as a tool to evaluate the capacity for a stormwater pond as a candidate for reuse and/or the need for the modification/retrofit. In addition, it may provide designers and planners with guidance to define standards for stormwater reuse.

#### **2. Materials**

#### *2.1. Study Area*

The City of Calgary is a semi-arid city in Alberta, Canada, with many stormwater ponds that are being considered as candidate sources of stormwater reuse for irrigating parkland and other public lands during the irrigation season [20]. The Inverness pond is one of the largest stormwater ponds in Calgary and is located in the southeast quadrant of the city. It contains approximately 235,000 m3 of water at the permanent water level. Seven inlets convey stormwater runoff from 415 ha of catchment area into the pond. Figure 1 shows the pond and the locations of the inlets and outlets. All of the inlets are submerged except for I1. Outlet O1 at the west wing of the pond is the main outlet of the pond and is a 1.5 m diameter concrete pipe. I4 inlet, located at the south wing, conveys water from the largest subbasin of the catchment, which is 257.97 ha. The second largest subbasin corresponds to I3 inlet with an area of 89.52 ha. In addition, there are four sediment forebays located in front of I2, I3, I4 and I5 inlets. Table 1 shows the subbasin area and design parameters corresponding to each inlet.


**Table 1.** Inlet and outlet parameters of the Inverness stormwater pond.

<sup>1</sup> Permanent water level.

**Figure 1.** Arial schematic of Inverness pond showing the location of inlets and outlets (base map from Google Earth (50◦54 39.4" N 113◦57 46.2" W at white cross-hairs at center), arrows and texts added by author).

#### *2.2. Data Collection*

Over a three-year data collection campaign, several data were collected within and around the pond. Water quality and flow rate data collected during several rain events in the 2007 irrigation season were taken at a skimming manhole located just upstream of the I5 inlet. The water quality data included FIB level (*E. coli*, FC, and total coliform), total suspended solid (TSS), turbidity, pH, and dissolved oxygen (DO) [20]. Rainfall data in 5-min intervals were acquired from the City of Calgary (rain gauge #26 located 1 km north of the study site in McKenzie Towne) and wind data were extracted from the Calgary International Airport. Five rain events documented in 2007 were used in this research. Three of these were lower intensity events leading to moderate flow rates (on 28 May, 26 August and 20 September 2007), while the other two (6 June and 12 September 2007) were high-intensity events leading to large flow rates into the pond. The rain event data are tabulated in Table 2.



An autosampler with 24 pre-sterilized bottles was placed inside the skimming manhole and collected samples for bacteriological analysis once triggered by a rain gauge. The autosampler was programmed to collect 12 samples (two bottles per sample) during long events, and less than 12 samples were analyzed for short duration events. The samples were collected at unequal intervals (3 to 50 min intervals), and much attention was paid to collect samples more frequently in the early stages of the rain event in order to catch the first flush effect. The water samples were recovered right after rain events or very early the next morning in case rain events occurred at night. The samples were packed with ice and analyzed at the Alberta Provincial Laboratory for Public Health by assaying microorganisms using the membrane method [20].

#### **3. Methods**

The model is comprised of a hydrological component that simulates stormwater runoff feed into the CFD, which simulates the hydrodynamics and predicts the bacteria levels within a water body. The modeling approach was formulated in a way in which *E. coli* are primarily discharged into the pond through contaminated stormwater runoff (water fowl droppings from the surface were not considered) and *E. coli* concentrations were affected by the biological fate of *E. coli* and pond hydrodynamics, which were in turn affected by meteorology and pond bathymetry. Studies in the literature have revealed that the concentration of bacteria is highly dependent on land use-in particular, the percentage of imperviousness [21]. Developed areas, which have more impervious areas, had elevated increases in *E. coli* concentrations in runoff [22] as opposed to less developed areas. In order to provide insight into the relationship between land use and bacteria concentration in runoff, bacteria concentrations were measured during storms in different watersheds with different land uses [23]. They demonstrated the variability in fecal coliform concentrations for different land uses. In the present study, the fraction of bacteria in stormflow corresponding to each land use was calculated using the findings of Schoonover and Lockaby [23]. However, here, all of the land uses except for newly graded land and farmland were considered as urban. Thus, by knowing the fraction of bacteria that could be generated from the area draining into inlet I5, the bacteria load of different inlets was estimated based on the measured data of inlet I5.

#### *3.1. Hydrological Model, Calibration and Validation*

The Hydrologic Modeling System (HEC-HMS 4.2.1) (Hydrologic Engineering Center, Davis, CA, USA) [24] was used to simulate stormwater runoff generated by each subbasin and drained into the pond in storm events. In the HEC-HMS, the Soil Conservation Service (SCS) curve number (CN) loss method [25], which is a very common and simple method for runoff estimation [26] was adopted. The SCS-CN model considers precipitation excess as a function of cumulative precipitation, soil cover, land use, and antecedent moisture. Precipitation excess was estimated by Equation (1):

$$P\_{\varepsilon} = \frac{\left(P - I\_{a}\right)^{2}}{P - I\_{a} + S} \tag{1}$$

where *Pe* is the accumulated rainfall excess (mm) at time *t*; *P* is accumulated rainfall depth (mm) at time *t*; *S* is the potential maximum retention (mm); *Ia* is the initial abstraction (mm), where *Ia* = 0.2*S*. Here, *S* is a function of the *CN*, which is a function of watershed characteristics, and was calculated by Equation (2) [27]:

$$S = \frac{25 \text{A}00 - 25 \text{C} \text{N}}{\text{CN}} \tag{2}$$

The *CN* number of each subbasin was initially estimated based on the land uses within each subbasin. The time step of the HEC-HMS model was set at 5 min. The runoff at I5, measured in the event on 6 June 2007, was employed to adjust the initial *CN* of this subbasin aiming to reproduce the runoff well. Accordingly, the *CN* of other subbasins was modified using the difference between the adjusted *CN* number for I5 subbasin and its initial *CN*.

#### *3.2. CFD Simulations*

This paper adopted Reynolds averaged Navier Stokes (RANS) equations as follows:

$$\rho \left[ \frac{\partial \overline{u}\_i}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\overline{u}\_i \overline{u}\_j) \right] = -\frac{\partial p}{\partial \mathbf{x}\_i} + \rho \frac{\partial}{\partial \mathbf{x}\_j} \left( \nu \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} \right) - \frac{\partial}{\partial \mathbf{x}\_j} \left( \rho \overline{u'\_i u'\_j} \right) \tag{3}$$

$$
\overline{u} = \frac{1}{T} \int\_0^T udt\tag{4}
$$

where the velocity decomposed into mean velocity and velocity fluctuation (*ui* = *ui* + *u i* ); *T* is period of time (*s*); ρ is density (kg/m3); *t* is time; *p* is pressure (pa); υ is kinematic viscosity (m2/s); *xj* is an axis of the Cartesian coordinate system, and *i,j* = 1,2,3 are indices indicating the three axes in the coordinate system. The term −ρ*u i u <sup>j</sup>* is called Reynold's stress term [28]. Since the popular *k-*ε turbulence model has been successfully applied to simulate water body flow fields [29,30] in this study the *k-*ε turbulence model [31] was used to calculate the Reynold's stress terms. The turbulence model is not detailed here. The fundamentals of CFD and turbulence modeling can be found in Versteeg and Malalasekera [32].

Kunkel et al. [33] studied the attachment of *E. coli* to particles in stormwater. The study showed that *E. coli* appeared to attach predominantly to fine particles (<4 μm), while a further study on attachment to smaller particles was recommended [33]. Another study on the attachment of *E. coli* showed that most of them attached to particles smaller than 2 μm [34]. The City of Calgary characterizes the particles in stormwater and suggests the settling velocity of 0.00592 (mm/s) for small particles (<10 μm) [1]. Using Stokes relation for settling velocity [35], the settling velocity of finer particles (<2 μm) can be estimated using:

$$w\_{\circ} = \frac{\mathcal{g}(\rho\_p - \rho\_w)}{18\rho\_{w\circ}\nu} d\_p^2 \tag{5}$$

where *ws* is settling velocity (m/s); ρ*<sup>p</sup>* and ρ*<sup>f</sup>* are particle and fluid density (kg/m3), respectively; *dp* is particle diameter (m); *g* is gravitational acceleration (m/s2). This relation is valid for small Reynolds numbers (*Re* < 1) [36]. With regard to the settling velocity of small particles, this condition was certainly satisfied in the Inverness pond [1]. The calculated settling velocity for the very fine particles (<2 μm) was 0.001 mm/s using the Stokes relation. This suggests that the attached bacteria settle less than 9 cm in 24 h. Therefore, it was assumed that all the bacteria (i.e., free-floating and attached *E. coli*) remain in the water column during the simulation. Accordingly, bacteria (both attached and free-floating bacteria) were modeled as passive scalars, which are massless particles that do not affect the physical properties of the flow field.

The passive scalar transport equation was solved for the transport of both free-floating and attached *E. coli*. The transport equation for a passive scalar, φ*<sup>i</sup>* is:

$$\frac{\partial}{\partial t} \int\_{V} \rho \phi\_{l} dV + \oint\_{A} \rho \phi\_{l} \upsilon. dA = \oint\_{A} j\_{l} dA + \int\_{V} S\_{\phi\_{l}} dV \tag{6}$$

where *t* is time (s); *V* is volume (m3); *A* is area (m2); *i* is the component index; *ji* is diffusion flux (kg/m2·s); *<sup>S</sup>*φ*<sup>i</sup>* is source term (kg/m2·s).

The fate of bacteria has been modeled using a negative value for the source term in the passive scalar transport equation in previous studies on waste stabilization ponds [18,19]. This approach is, however, not appropriate for stormwater ponds due to two reasons. Firstly, the flow field and bacteria transport in stormwater ponds are intrinsically unsteady; thus, the location and strength of the sources cannot be determined. Secondly, bacteria concentrations in some spots of stormwater ponds are very low, so a general source, representing decay over the entire domain, might cause negative bacteria concentrations which are not physically feasible. Therefore, in this paper, the transport of bacteria was calculated with Equation (6), while the fate of bacteria was modeled separately.

A finite volume discretization scheme was used to discretize Equations (3) and (6), and the turbulence model equations over the domain to solve the flow field and the concentration of passive scalar. Then, a field function was implemented for the first-order kinetic decay (Equation (7)) for the free-floating portion of bacteria and added the concentrations of both attached and the remaining free-floating bacteria at each computational cell. Characklis et al. [13] found that 30–55% of bacteria attach to sediments in stormwater, whereas the other study [37] applied a higher value in the range (50%). Therefore, in this paper, it was assumed that 50% of bacteria are attached to sediments and do not decay. Decay is calculated using Equation (7).

$$N = N\_0 \exp(-\mu t) \tag{7}$$

where *N* and *N*<sup>0</sup> are the numbers of indicator bacteria at time *t* and at *t* = 0, respectively; μ is the decay rate (*s*−1) [38]. The decay rate of *E. coli* was calculated with the following relationship that has been used in several studies [14,39] in the literature:

$$\mu = k\_{20} \frac{e^{\left(-\frac{(T-25)^2}{400}\right)}}{e^{\left(-\frac{25}{400}\right)}}\tag{8}$$

where *<sup>k</sup>*<sup>20</sup> is the decay rate of *E. coli* at 20 ◦C and *<sup>T</sup>* is water temperature (◦C). A value of 1.25 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>s</sup>−<sup>1</sup> was proposed for *k*<sup>20</sup> of *E. coli* in the water column; average measured data were used for temperature in this paper.

#### 3.2.1. Boundary and Initial Conditions

The effect of wind on the flow field was not negligible considering the size of the study pond. Since the wind vector was not stationary with time, a dynamic boundary condition of wind on the pond surface was defined. The effect of the wind was applied as a velocity vector on the top layer of the pond. The velocity of the top layer of the pond was calculated from the wind velocity based on the experimental data from a study by Banner and Peirson [40], which measured velocity at the top layers of water bodies under windy conditions.

The boundary condition of the inlets was set as the "velocity inlet" with parameters changing with time including velocity and *E. coli* concentration. A time variable boundary condition was also defined which linearly interpolated the value of bacteria concentration and velocity between the interval that data were collected. "Pressure outlet" and "stationary wall" were set as the boundary conditions for the outlets and the bottom of the pond, respectively.

A steady-state simulation in no storm conditions and with an average wind was simulated and its result (flow field) was used as the initial condition for the main simulations. The main simulation is unsteady. Using the steady simulation's results also accelerated the convergence. The unsteady simulations started one hour before the storm in order to let the flow field form based on the real wind data before the storm.

#### 3.2.2. Other CFD Settings

The domain of the CFD simulation was sketched based on the bathymetric data of the pond comprised of approximately 33000 GIS points. AutoCAD CIVIL 3D 2018 was used for sketching the domain that includes the water body of the pond and a few meters (10 m on average) of inlet and outlet pipes. A grid was made that is more highly clustered close to the inlets, outlets and water edges. In the CFD simulations, grid dependence was checked by comparing velocity from the simulation with different numbers of grids. An unstructured grid with hexahedral cells and a number of 1.5 million was found sufficiently fine and the corresponding results were presented. Using finite volume discretization, RANS equations with the *k-*ε turbulence model were solved in an unsteady-state with the time step of 0.5 s and 24 inner iterations between each time step. The time step was selected by an approach

similar to the grid independence comparison but the number of inner iterations was chosen based on the convergence trend of the solution. The pond was simulated from 1 h before the start of an event until 24 h after the end of the event. The simulations were run on Cedar High-Performance Computer, located at Simon Fraser University, Vancouver, BC, Canada, with 576 GB of memory and 192 computational cores (each core has two CPUs of 2.1 GHz). The CFD simulations were performed by using commercial CFD code STAR-CCM+ 12.04.011 [41].

#### **4. Results and Discussions**

#### *4.1. Hydrological Modeling Performance*

The hydrological model was manually calibrated with the measured data in inlet I5 on 6 June 2007 and then validated on four other events: 28 May, 26 August, 12 September, and 20 September 2007. During the calibration, the error in the total volume of water between prediction and observation was adopted to evaluate model performance, as the total volume of water entering the pond was considered to be the key factor that determines bacteria loading.

The CN number for the subbasin draining to I5 was calculated to be 11% less than the initially calculated CN value in the calibration of the 6 June 2007 event. Noting this, all other CN numbers were decreased by 11%. Table 3 shows the CN numbers for the different subbasins before and after calibration of the model. Lag time was also calibrated at the same time, and a 30 min lag time was found to be optimal. In addition, an attempt was made to calibrate initial abstraction; however, the default setting (*Ia* = 0.2*S)* resulted in the lowest error and was, therefore, left unchanged. Figures 2–5 show the modeled and observed hydrographs for the validation events for subbasin I5. The Nash-Sutcliffe model efficiency (NSE), which is a reliable criterion for assessing the goodness of fit of hydrological models [42], and the error in total volume for all events, are tabulated in Table 4. The calibrated CNs (Table 3) were then used to model other subbasins.

**Table 3.** Curve numbers of the different subbasins draining to the Inverness pond inlets.


**Figure 2.** Modeled and observed hydrographs at I5 inlet of storm event on 28 May 2007.

**Figure 3.** Modeled and observed hydrographs at I5 inlet of storm event on 26 August 2007.

**Figure 4.** Modeled and observed hydrographs at I5 inlet of storm event on 12 September 2007.

**Figure 5.** Modeled and observed hydrographs at I5 inlet of storm event on 20 September 2007.


**Table 4.** Model performance in model calibration and validation. The event on 6 June 2007 is the calibration event and all others are validation events.

Results show good agreement between the modeled and measured data. In all of the events except for 20 September 2007, the model slightly overestimated the total volume but the total volume error is less than 10%. The underestimation of the total volume for the event on 20 September 2007 may be due to the excessive moisture in the soil resulting from heavy rains that fell prior to the event; thus, potentially leading to greater runoff generation. However, the Nash-Sutcliffe efficiency in all of the modeled events is reasonable.

#### *4.2. CFD Simulations*

The simulation results for the event on 26 August 2007 are detailed and discussed in this section as an example. The concentrations of *E. coli* on the surface of the pond at different time steps after the rain event are shown in Figure 6. As time passes, bacteria entering from the inlets redistribute in the pond. The flow field was affected by the inlet velocities for the first few hours after the events, but afterwards, the wind was the only parameter affecting the flow field.

**Figure 6.** Contours of *E. coli* concentration (cfu/100 mL) on the surface of the Inverness pond on (**a**) 26 August 2007 at 5 p.m., (**b**) 26 August 2007 at 11 p.m., (**c**) 27 August 2007 at 5 a.m., and (**d**) 27 August 2007 at 11 a.m.

The bacteria concentrations in the south wing of the pond were greater as compared with the other two wings (the east and west wings). There are two inlets (I5 and I4) in the south wing. The *E. coli* concentration is higher in the stormwater runoff entering the pond from inlet I5 than all other inlets; while bacteria loading from I4 is highest as it drains the largest area (Table 1). For example, the *E. coli* concentration in the storm runoff from I5 on 26 August 2007 at 2 p.m. (during the storm) was 2100 cfu/100mL, while the storm runoff from I4 had a concentration of 1 038 cfu/100 mL. However, the flow rates of the inlets at that time were 0.045 m3/s and 0.39 m3/s for the I5 and I4 inlets, respectively. Therefore, more bacteria mass entered the pond from inlet I4 even though the concentration of bacteria was greater in I5. In addition, a relatively large amount of bacteria entered the pond from inlets I3, I7 and I6, but most of the bacteria that came from I7 and I6 immediately exited the pond through outlet O1, because these two inlets are in proximity to the outlet. Therefore, they do not have a significant effect on the bacteria level in the pond. In general, inlets I4, I3 and I5 had the most significant effect on the bacteria levels in the pond.

Figure 7 shows the vertical profile of *E. coli* distribution in the pond at 6 h after the end of the storm. It reveals that the bacteria concentration also changes with depth. As illustrated in Figure 7, the maximum concentration of *E. coli* on the surface barely reached 90 cfu/100 mL. However, the maximum *E. coli* concentration was more than 120 cfu/100 mL at 2 m depth. In addition, the bacteria in the middle of the pond (where the three wings join) was less distributed at the bottom compared to the surface. The reason is that the direction of the wind had been NE and NNE for the last few hours, so the bacteria escaping from the sediment forebay of inlet I2 could reach the other side where the south wing and the west wing join. In contrast, there was current flow toward the NE and NNE directions at the bottom simply because of mass conservation. This current potentially brought clean water close to the I2 sediment forebay. There was also a downwelling near the bank of the middle of the pond where the south and the west wings join. Generally, the wind causes differences in *E. coli* distribution in different layers of the water column. *E. coli* data collected at different depths on five random days between the year 2006 and 2007 [2] reveal that the bacteria concentration changed with depth. However, these data did not show any specific trend, and this is likely an indication of the influence of multiple environmental conditions on the bacteria distribution at various depths.

**Figure 7.** The vertical profile of *E. coli* concentration on 26 August 2007 at 11 p.m. (**a**) on the surface, (**b**) at 0.5 m below the surface, (**c**) at 1 m below the surface, and (**d**) at 2 m below the surface.

The sediment forebay corresponding to I3 was full of bacteria at all depths after the event (Figure 7). The sediment forebay was able to retain bacteria many hours after the event, which resulted in keeping the east wing relatively clean as compared with the south wing. Figure 8 magnifies the contour of *E. coli* concentration near the sediment forebays at the surface on 27 August 2007 at 2 a.m. It also confirms that the sediment forebay of I3 outperforms the forebay of other inlets. As illustrated in Figure 8, the

concentration of *E. coli* in the sediment forebays of inlets I4 and I5 appears to be similar or even slightly less than that in the region outside the forebays. This reveals that the bacteria entering these forebays during the storm were promptly discharged out of the forebays and consequently increased the *E. coli* concentration in their nearby regions.

**Figure 8.** *E. coli* distribution on the surface 27 August 2007 at 2 a.m. (9 h after the end of the storm).

The design of sediment forebays, such as their configuration and size, determines their efficiency in trapping bacteria. Figure 9 shows streamlines coming out from the inlets and spreading throughout the pond during the storm. The streamlines from I4 continue straight out of the forebay, which means the bacteria were transported into the pond. One reason is that the size of the forebay was not large enough to fit a circulation proportional to the high flow rate of the I4 inlet. Another reason may be the configuration of the inlet and sediment forebay. The direction of the streamline from I4 and the corresponding sediment forebay were perpendicular to each other (i.e., they were in front of each other), so the water jet coming out of the inlet easily escaped the forebay without circulating. The situation was the same in the sediment forebays corresponding to the inlets I2 and I5. However, their size was proportional to their flow rate. The sediment forebay corresponding to the inlet I3 had the best configuration since the direction of streamlines coming out of the inlet was parallel to the forebay. Thus, the bacteria coming from the inlet had no way but to circulate and in the long term, they die, which kept the east wing relatively clean. In addition, the size of the forebay was large enough to fit two large eddies.

During the data collection campaign, surface grab samples were collected from six different sites at the pond over several days [2]. All of the samples were collected at an average depth of 15 cm. One of the sample collection days was after a day with heavy rain. On 10 September 2005, 68 mm of rain fell and on the next day, water samples were collected from the six sites. Although much greater rain fell on 10 September 2005 than in the simulated events from 2007, it provided an interesting case for validation due to the dominant effect of rain on the bacteria distribution. The initial bacteria level in the pond before an event can influence the final bacteria levels after an event, particularly for small, low rainfall events. However, for high rainfall events, the majority of bacteria is transported into the pond with the storm runoff, and the initial bacteria level is a far smaller proportion of the total bacteria level. In this paper, it was assumed that the initial bacteria level was negligible, and only the effect of the storm was studied. Thus, the modelling results can be compared with the collected data on 10 September 2005. Although it was noted that due to a lack of data it was difficult to validate the process-based FIB methods thoroughly [11], a comparison was done between the different events using a non-dimensional number. The non-dimensional number was computed as the ratio of *E. coli*

concentration at a site to the maximum *E. coli* concentration among all of the six sites at a certain time. Thus, the non-dimensional number at the site where *E. coli* concentration was a maximum was one. The number was calculated for the individual modeled events, and then the average was taken for each site. Figure 10 shows the non-dimensional number at the six sites for the measured data and the average of the simulated events with its variation range. Simulation results show a good agreement with the measured data in recognizing hot spots (spots with a high concentration of bacteria) and spots with the lowest level of bacteria. The variation in the non-dimensional number in the east and west wings were higher than that in the south wing and in the middle of the pond. This may be due to a stronger influence from meteorological factors such as wind and rain at the tip of the east wing and the tip of the west wing. On the other hand, the low variation in the calculated data in the south wing and their high values show that the south wing always had the highest bacteria levels in the pond. Therefore, it is not recommended to extract water for reuse from the south wing.

**Figure 9.** Velocity streamlines on 26 August 2007 at 3 p.m. (during the storm).

In addition, both measured and calculated data show that the middle of the pond, where the three wings join, had the lowest level of bacteria among the six sites; this site was a reliable spot for water extraction due to low concentrations of bacteria and low variability in the calculated data. There were no collected data for the middle of the west wing; however, simulation results show that this spot had the lowest velocity and *E. coli* concentrations in the pond. In a region with low velocity, fewer bacteria can enter. In addition, even if some bacteria enter a low-velocity region, most of them die due to high retention time in the region, which keep the bacteria concentration relatively low. In addition, the west wing has two inlets and an outlet all located at the tip of the wing. Thus, most of the bacteria that were discharged to the wing by the inlets exited the pond promptly. Therefore, the middle of the west wing can be considered as a potential location for water extraction. However, more data should be collected from the entire pond for a more comprehensive validation. In addition, further modeling studies after rain events are recommended when more data are available. More experimental and modeling studies need to be conducted in order to study the effect of sedimentation on bacteria transport in the pond, which was assumed to be negligible in this paper.

**Figure 10.** Non-dimensional *E. coli* concentration 15 cm below the surface.

#### **5. Conclusions**

An integrated hydrological-CFD model was developed to simulate bacteria concentrations within the Inverness stormwater pond in Calgary, Alberta, in order to locate the best (cleanest) locations within the pond to extract water for potential reuse. The integrated model was calibrated at one inlet with observed flow rate data in order to estimate the flow rates and bacteria levels at other inlets. The flow rate and bacteria concentration at the inlets were implemented as boundary conditions for an unsteady CFD simulation to determine the distribution of bacteria in the pond. The bacteria concentration distribution in Inverness pond was simulated both during and after four rain events using the integrated model. Results showed a good agreement with collected data from six different points in the pond, indicating that the model can be successfully used to model fate and transport of bacteria in stormwater ponds. It was found that the wind plays a crucial role in forming the flow field in the pond, which affects the bacteria distribution. In addition, one of the sediment forebays, located in the east wing of the pond, successfully trapped large amounts of bacteria until death occurred. However, other sediment forebays were found ineffective. The middle of the west wing was found to be the best location for extracting water because it showed the lowest levels of bacteria during the simulation of the pond both in and after the events.

**Author Contributions:** Conceptualization, modeling, simulation, writing, F.A.; supervision, investigation, reviewing, editing, and fund acquisition, C.V., J.H.; fund acquisition, project administration, and reviewing, N.F.N.

**Funding:** This research was funded by Alberta Innovates, grant number RES#0030866.

**Acknowledgments:** The authors would like to thank Mostafa Rahimpour for his comments on CFD modeling and Belaid Moa from WestGrid for his technical support on the parallel computation of the work.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Near-Wake Flow Structure of a Suspended Cylindrical Canopy Patch**

#### **Ay¸se Yüksel Ozan 1,\* and Didem Yılmazer <sup>2</sup>**


Received: 26 November 2019; Accepted: 21 December 2019; Published: 25 December 2019

**Abstract:** Urban stormwater is an important environmental problem, especially for metropolitans worldwide. The most important issue behind this problem is the need to find green infrastructure solutions, which provide water treatment and retention. Floating treatment wetlands, which are porous patches that continue down from the free-surface with a gap between the patch and bed, are innovative instruments for nutrient management in lakes, ponds, and slow-flowing waters. Suspended cylindrical vegetation patches in open channels affect the flow dramatically, which causes a deviation from the logarithmic law. This study considered the velocity measurements along the flow depth, at the axis of the patch, and at the near-wake region of the canopy, for different submerged ratios with different patch porosities. The results of this experimental study provide a comprehensive picture of the effects of different submergence ratios and different porosities on the flow field at the near-wake region of the suspended vegetation patch. The flow field was described with velocity and turbulence distributions along the axis of the patch, both upstream and downstream of the vegetation patch. Mainly, it was found that suspended porous canopy patches with a certain range of densities (SVF20 and SVF36 corresponded to a high density of patches in this study) have considerable impacts on the flow structure, and to a lesser extent, individual patch elements also have a crucial role.

**Keywords:** suspended vegetation; FTW; ADV; velocity profile; submerge ratio; SVF

#### **1. Introduction**

Increases in urban and agricultural development cause stormwater runoffs that contains nutrients that can have severe effects on aquatic ecosystems and human health [1]. The use of floating treatment wetlands (FTWs) is a new and innovative tool to remove nutrients in stormwater wet detention ponds [2]. FTWs represent a human-made ecosystem similar to natural wetlands [3]. A typical FTW is shown in Figure 1 [4].

**Figure 1.** An example of floating treatment wetland [4].

Vegetation patches or layers are essential in aquatic environments in submerged, emergent, or suspended forms. Suspended canopies, which are the focus of this study, form on the water surface, with stems toward the bottom that do not touch the bottom, in contrast to submerged and emerged canopies that are attached to the bottom. *Lemnaceae, Eich- hornia crassipes, Pistia stratiotes, Salvinia molesta,* and *Laminariales* are typical examples of suspended canopies [5–7]. *Lemnaceae* are used extensively for the purpose of purifying wastewater and eutrophic water bodies [8]. Suspended canopies in flow have been investigated in terms of the ecological effects [9–12] and the mean velocity and turbulence structures of open-channel flows included in these types of canopies [6,7]. Fishing activities and marine navigation can be at risk because of these types of floating vegetation patches [8]. Bed friction has a substantial effect on the velocity and turbulence below and within floating vegetation patches [6,13], however, it has a small role on the flow structure inside submerged or emergent canopies [14–21]. Despite the increasing importance of suspended canopies, there is still a limited number of studies on them. Detailed flow structures in the wake regions of the suspended canopies should be examined to understand the transport of fluid and pollutants.

Liu et al. [22] studied longitudinal dispersion in flow along the floating vegetation patch. They also indicated that a large number of studies on the longitudinal dispersion in open-channel flows with submerged and emergent vegetation have been presented [16,23–25]. However, there is an insufficient number of studies on the effects of floating vegetation patches in the literature [22]. They proposed a four-zone model with which they defined the velocity profile and turbulent diffusion coefficient profile for the four zones in their study.

Huai et al. [7] studied the effect of suspended vegetation on velocity distribution at open-channel flows. They considered four different submergence ratios (hv/H = 0.125–0.5, where hv is submergence ratio of the suspended vegetation and H is the flow depth) with a constant density (defined as the canopy projected area per unit volume, a = 1.272 m<sup>−</sup>1) of the canopy in their experiments. They defined four parts along the flow depth with the suspended vegetation (Figure 2) and determined the velocity profile for every part with the solution of the momentum equation. They called Zone 1 and Zone 2 the outer and internal vegetation layers, respectively. They explained that the drag force of vegetation is dominant in Zone 1 where the velocity is slow and almost constant. The velocity in this zone corresponds to that of the wake region in the submerged vegetated flow [7]. They defined the non-vegetated layer with two zones: Zone 3 and Zone 4. Kelvin–Helmholtz (K–H) vortices characterize the middle mixing layer (tml), which is placed near the bottom of the suspended vegetation [22]. The K–H vortices control the vertical mass transport in the lower region of the suspended vegetation and in the upper part of the gap, which is the non-vegetated region of the flow [26,27]. Zone 2 corresponds to the penetration depth, hi, where the length of the K–H vortices penetrate the vegetation region [22]. The non-vegetated area in the mixing layer is defined as Zone 3, in which the thickness is *tml*−*hi* [22]. The resistance caused by suspended vegetation affects Zone 1, Zone 2, and Zone 3, while the friction resistance of the bed affects Zone 4, as it is located close to the bed [22]. The velocity has a maximum value at the junction of Zones 3 and 4. Liu et.al. [22] stated that because of the maximum velocity, the Reynolds shear stress has zero value at this point where the impacts of river bed and suspended vegetation are consistent.

**Figure 2.** Schematic velocity distribution in suspended vegetated flow [22].

The wake structure downstream of a circular suspended porous obstruction, which represents a model of a suspended vegetation patch, was investigated in this paper. The main question of this study is how the submergence ratio affects the structure of the wake and flow close to the bed, downstream of the suspended patch. For this purpose, the velocity measurements were performed along with the flow depth downstream of the canopy for different submergence ratios with different patch porosities.

#### **2. Experimental Methods**

Experiments were performed in an 11 m long, 1.2 m wide, and 0.75 m deep re-circulating glass flume with a concrete horizontal bed at the Civil Engineering Laboratory of Aydın Adnan Menderes University (Figure 3). The suspended vegetation patch was of diameter D = 0.3 m and was formed by rigid plastic cylinders of diameter d = 0.01 m (Figure 4). The cylinders were arrayed in a staggered pattern (Figure 5). The suspended patch was placed at the center of the flume in the transverse direction. The number of cylinders per unit bed area, n (IP/cm2), defines the density of the patch. Here, IP means individual plant number per square centimeter and area is the A = πD2/4. Vegetation density is defined as the frontal area per unit volume, given with a = nd (cm<sup>−</sup>1), and characterizes the blockage caused by the vegetation layer [28]. The average solid volume fraction is defined with <sup>φ</sup> = nπd2/4 <sup>≈</sup> ad [29]. In the experiments, three different densities of the canopy patch and a solid cylinder were used to simulate the suspended vegetation (Table 1). The Plexiglas solid cylinder was composed as a closed cylinder with D = 30 cm to represent the suspended Solid Case (SC) in the flow area (Table 1). Additionally, a constant mean depth-averaged velocity (*Udm* = 0.116 m/s) and constant flow depth (H = 0.3 m) were considered in the experiments. Here, velocity values were integrated along the flow depth to calculate mean depth-averaged velocity (*Udm*). Mean depth-averaged velocity (*Udm*) was used to calculate Reynolds number (Re = *UdmH*/ν) and Froude Number *Fr* = *Udm*/ *gH* for the non-vegetated case. Here, *g* is the gravitational acceleration, and ν is the kinematic viscosity. Patches were immersed in the flow at two different heights: hv1 = 0.1 m and hv2 = 0.2 m. The cases at the experiments are outlined in Table 1.

**Figure 3.** Experimental set-up: (**a**) plan view, (**b**) cross-section (units are in cm).

**Figure 4.** *Cont*.

**Figure 4.** Photos from vegetation patches for different SVFs (Solid Volume Fractions): (**a**) side view, (**b**) top view, (**c**) view from an experiment with an Acoustic Doppler Velocimeter (ADV).

**Figure 5.** Schematic sketch of the plants of different densities; (**a**) φ = 0.05, (**b**) φ = 0.20, (**c**) φ = 0.36, (**d**) Solid Case (units are in cm).


**Table 1.** Experimental conditions.

In this study, *x*, *y*, and *z* are the longitudinal, transverse, and vertical coordinates, respectively. The origin of the coordinate system is at the upstream edge of the suspended vegetation for the *x–y* plane and the flume bed for the *z* direction (Figure 6). A 3-D Acoustic Doppler Velocimeter (ADV) SonTek 10 MHz with a 3-D down-looking probe was used to measure three instantaneous components of velocity (u, v, w), corresponding to *x*, *y*, and *z* directions shown in Figure 6, respectively, during 120 s at 25 Hz. The velocity range was ±0.30 m/s, with a measured velocity accuracy of ±0.5%, and the sampling volume diameter was 6 mm established 0.05 m down from the probe in the experiments. The velocity measurements were conducted at 11 points in the x-direction and at eight different levels along the flow depth (Figure 6). Measurement points along the *x*-direction were chosen by moving at intervals of 15 or 30 cm upstream and downstream of the vegetation patch (Figure 6c). Alongside this, measurement points along the Z direction were placed at different intervals which correspond to

Z1 = 0.5 cm, Z2 = 1.0 cm, Z3 = 1.5 cm, Z4 = 2.0 cm, Z5 = 3.0 cm, Z6 = 4.5 cm, Z7 = 6.0 cm, Z8 = 8.0 cm, Z9 = 11.0 cm, and Z <sup>10</sup> = 15.0 cm (Figure 6a,b). The ADV company suggests that a 15 db signal-to-noise ratio (SNR) and a correlation coefficient larger than 70% for high-resolution measurements [30]. For this reason, conditions where SNR <15 db and the correlation coefficient <70% are used to remove the low-quality samples in the measured data.

**Figure 6.** Velocity measurement points (**a**) along the flow depth for hv/H = 0.33, (**b**) along the flow depth for hv/H = 0.67, and (**c**) along the flow direction (units are in cm).

#### **3. Results and Discussion**

The longitudinal profiles of streamwise averaged velocity, U, for various different densities of patches at mid-depth (z/H = 0.5), for hv/H = 0.67, and hv/H = 0.33, are presented in Figure 7. Here, the vegetated part (VP) corresponds to the suspended vegetation canopy and the mid-depth corresponds to beneath the patch for hv/H = 0.33. The diversion of flow begins almost one diameter upstream of the patch, which is consistent with Rominger and Nepf's [31] results (Figure 7a). There is a recirculation zone (U < 0) downstream of the flow, which extends to 1.5 D and 2 D distance from the edge of the patch for SVF20 and SVF36, respectively (Figure 7a). As the density increases, the recirculation zone moves downstream and the magnitude of negative streamwise velocity increases. Moreover, there is no recirculation zone for the lowest density (SVF5). Zong and Nepf [29] found that there is a recirculation region in the wakes behind denser patches (φ = 0.10), but absent in the wakes of sparser patches (φ = 0.03). This result is consistent with our results. A steady wake region was observed downstream of the patch, close to the patch for SVF5, which is not present at higher SVFs in this study. Here, steady wake corresponds to 2 D distance beyond the end of the patch, where U velocity is almost constant.

**Figure 7.** Longitudinal profiles of U velocity for different densities of patches at mid-depth (z/H = 0.5) for (**a**) h/H = 0.67 and (**b**) h/H = 0.33.

Zong and Nepf [29] mentioned that the wake can never perfectly recover since the drag-induced momentum deficit should be constant downstream of the obstruction. They further pointed out that the total length of the wake, L, is defined as the distance from the downstream end of the patch to a point where the velocity recovery rate is induced to (∂(*U*/*Udm*)/∂(*x*/*D*) < 0.1). The total length of the wake (L) is almost 3 D at mid-depth for the solid body and for ReD <sup>=</sup> 3.5 <sup>×</sup> 104 (Figure 7a). This value is consistent with Zong and Nepf's [29] results, in which L = 3 D and L = 3.3 D for ReD = 2.2 <sup>×</sup> 104 and ReD = 4.1 <sup>×</sup> 104, respectively, in their study. Cantwell and Coles [32] also defined L = 2.5 D for a cylinder wake at ReD <sup>=</sup> 1.4 <sup>×</sup> 105 in their study. The recovery region extends to 5 D beyond the end of the patch in our study. Another noteworthy situation is that U increases rapidly to *Udm* at the highest density (SVF36) compared to the lower densities (SVF5 and SVF20) (Figure 7a).

Considering the lower depth ratio (h/H = 0.33), there is almost no effect of the patch beneath the patch at mid-depth for the solid body and lowest SVF (Figure 7b). However, suspended vegetation still affects the flow beneath the patch for a higher density of patches (SVF20 and SVF36). The obstruction effect on the flow is higher at the highest patch density (SVF36).

Longitudinal profiles of time-averaged transversal velocity, V, for different densities of patches and depth ratios at mid-depth were also examined (Figure 8). A suspended vegetation patch does not have a significant effect on transversal velocity beneath the patch for a low solid volume fraction (SVF5) at a lower depth ratio (hv/H = 0.33) (Figure 8a). As the suspended canopy become denser, the magnitude of the transversal velocity increases downstream until 2 D distance from the end of the patch for a lower depth ratio (hv/H = 0.33). Similar to a lower depth ratio, there is almost no effect on transversal velocity downstream of the patch for low solid volume fractions at higher depth ratios (hv/H = 0.67) (Figure 8a). There is diverging flow (*V* > 0) downstream of the patch which extends until 5 D distance for SVF20, 3 D distance for SVF36, and 1 D distance for the Solid Case. The position of diverging flow occurs closer to the patch with the increasing SVF. Transversal velocity varies around zero beyond the point where the diverging of the flow ends.

**Figure 8.** Longitudinal profiles of V velocity for different densities of patches and depth ratios at mid-depth (z/H = 0.5): (**a**) φ = 0.05, (**b**) φ = 0.20, (**c**) φ = 0.36, (**d**) Solid Case.

Figure 9 shows the vertical profiles of the normalized velocity U/*Udm* at different upstream and downstream positions (x'/D) for hv/H = 0.67. Here, x' refers to the distance from the patch edge and NV and NVlog correspond to the normalized velocity (U/*Udm*) distribution and normalized logarithmic velocity (Ulog/*Udm*) distribution, respectively, at non-vegetated flow conditions. You can find detailed information about the calculation of Ulog from Yılmazer et al. [17]. A combination of a channel bed and suspended patch causes maximum velocity at the gap, with a flow similar to a wall jet beneath the vegetation patch. Here, the layer between the bed and the maximum velocity point is called the inner layer, where wall shear is effective; the area above the point of the maximum velocity is called the outer layer, where patch shear is effective. Flow has a maximum velocity of almost 1.45–1.5 *Udm* beneath the patch and downstream for higher densities of patches (SVF20 and SVF36). Comparatively, the maximum velocity for the lowest density of patches (SVF5) is around 1.3 *Udm*, which is similar to the Solid Case. It can be concluded that there is a suspended patch effect that increases the velocity beneath the suspended patch for a certain range of SVF, which corresponds to high densities of patches. Another important point is the placing of the maximum velocity, which is z/H - 0.15 for SVF20 and SVF36, and z/H - 0.3 for SVF5 and the Solid Case. Inner layer thickness decreases with the increasing SVFs. The placement and the magnitude of the maximum velocity at the gap region shows similar behavior at SVF5 and the Solid Cases, which is similar for higher densities (SVF20 and SVF36). It can be concluded that the individual patch elements have an important role for lower densities, which is SVF5 in this study. U velocity has negative values downstream of the patch for the densities SVF20 and SVF36, and its value increases with increasing patch density at the near-wake region (~x'/D ≤

1.5 D). There is no reversed flow region for the lowest density patch (SVF5) and solid patch. Moreover, the flow behavior at the Solid Case is similar to the lower densities of the patch, e.g., SVF5. Maximum velocity value bigger than *Udm* is present along x' < 3 D, downstream of the gap region, which also corresponds to the distance where flow behavior is similar to wall jet flow.

**Figure 9.** Vertical profiles of U velocity for different densities of patches along streamwise direction for h/H = 0.67: (**a**) φ = 0.05, (**b**) φ = 0.20, (**c**) φ = 0.36, (**d**) Solid Case.

The vertical profiles of U velocity for different densities of patches are shown in Figure 10 for hv/H = 0.33. Velocity distribution and flow behavior deviate at x' = 0.5 D distance downstream of the patch, and is totally different in the case of a higher submergence ratio. The patch does not have a big effect on the flow structure along with the entire flow depth for this submergence ratio. The behavior of the velocity distribution for SVF20 and SVF36 is similar, as is the case for SVF5 and the solid patch. Flow has a maximum velocity of almost ~1.25 *Udm* for higher densities (SVF20 and SVF36) and around ~1.1 *Udm* for the lowest density (SVF5) and Solid Case. This result shows us that it is still possible to see the effect of density for high densities (SVF20 and SVF36 in this study) downstream of the gap region for lower submergence ratios. Downstream of the gap region, a maximum velocity value bigger than *Udm* is present along x' < 3 D for SVF36 and SC, and x' < 5 D for SVF5 and SVF20. It can be concluded

that a high density of the submerged vegetation patches still has an influence on the flow structure beneath the patch, even at low submergence ratios.

**Figure 10.** Vertical profiles of U velocity for different densities of patches along streamwise direction for h /H = 0.33: (**a**) φ = 0.05, (**b**) φ = 0.20, (**c**) φ = 0.36, (**d**) Solid Case.

Longitudinal profiles of *urms*/*Udm* and *vrms*/*Udm* for different densities of patches at different depths (z/H = 0.2 and z/H = 0.5) along the streamwise direction for hv/H = 0.67 are shown in Figure 11. Here, z/H = 0.2 corresponds to the gap region and z/H = 0.5 corresponds to the upstream and downstream parts of the patch. First of all, velocity fluctuations are higher upstream of the gap region than upstream of the patch for all cases. Moreover, there is a peak at x'= –0.5 D upstream of the gap region that is not present upstream of the patch. This shows us that the effect of a submerged obstruction starts in the upstream section, below the obstruction. Turbulence starts with a decaying trend downstream, and then it starts to increase for SVF5 and SVF20 cases at both levels. Specifically, the shear layer caused by the patch increases the turbulence level beneath the patch. Turbulence shows a decaying trend for SVF36, which is the highest density of patches, and the values of turbulence become closer at x' = 3 D distance for different depths (z/H = 0.2 and z/H = 0.5). This shows us that the individual stem turbulence is not significant for the highest patch density (SVF36), but is significant at lower patch densities (SVF5 and SVF20). Individual cylinders cause small-scale turbulence, which

results in a higher turbulence region behind the patch. The turbulence level has the highest value, especially for SVF20, in the near-wake region of the patch in this study. This is because of the significance of the individual cylinders for the patch intensity.

**Figure 11.** Longitudinal profiles of *urms*/*Udm* and *vrms*/*Udm* for different densities of patches along the streamwise direction for hv/H = 0.67: (**a**) φ = 0.05, (**b**) φ = 0.20, (**c**) φ = 0.36, (**d**) Solid Case.

Comparing the cases of SVF20 and SVF36, the turbulence level decreases with increasing patch density (SVF36) to a level close to that of the solid body at mid-depth. This shows that the suspended patch with the highest porosity patch (SVF36) behaves like a solid body. The turbulence level is higher downstream of the solid body than it is downstream of the gap region along x' = 2 D distance, beyond which point it is reversed (Figure 11d). After this point, the magnitude of the turbulence is similar for the highest density (SVF36) of patches and the solid body at mid-depth. However, the magnitude of the turbulence is higher for the solid body than the highest density (SVF36) of patches in the recovery region. This may occur because of the wake region behind the solid body. There is only one point, which is the closest point to the patch, (x = 0.5 D) downstream of the highest density patch (SVF36), where turbulence is higher than the solid body at mid-depth. This point corresponds to the wake region of the individual stem of the highest density patch. This shows us that the highest density patches behave like a solid body at the wake region, with the exception of very close to the patch (x ≤ 0.5D).

Figure 12a shows the spectral density variation of transverse velocity with frequency for Case SVF5 and the Solid Case. The eddy dominant frequency corresponds to the peak velocity spectra *Svvpeak*. Case SVF5 and the Solid Case are discussed as representative examples and for the comparison of the canopy and Solid Cases (Figure 12). The dimensionless Strouhal number (*St* = *fL*/*Udm*) defines the relationship between the characteristic length of a vortex and the frequency. Here, *f* is the vortex shedding frequency (dominant frequency), *L* is the characteristic vortex length, and *Udm* is the characteristic flow velocity. *St* is almost constant around *St* = 0.2 within a wide range of Reynolds numbers, between Re <sup>=</sup> 250 to 2 <sup>×</sup> 105 for flows passing a cylinder [33]. Zhao and Huai [34] assumed that this can also be used for the bundle of cylinders in their study. Therefore, we considered *St* = 0.2 for both Case SVF5 and the Solid Case. For Case SVF5 (x'/D = 0.5 and z/H = 0.5), the vortex shedding frequency is *f* = 0.08 Hz (Figure 12a) and the characteristic velocity *Udm* is approximately 0.116 m/s. The characteristic vortex length *L* is almost 0.293 m, which is almost equal to the diameter of the patch

(D = 0.3 m). For the Solid Case (x/D = 0.5 and z/H = 0.5), the vortex shedding frequency is *f* = 0.11 Hz (Figure 12a) and the characteristic velocity *Udm* is almost 0.116 m/s. The characteristic vortex length *L* is almost 0.211 m, which is smaller than the canopy diameter. When the magnitude of the peaks in Svv is compared, it is clear that the turbulence intensity of the wake-scale vortices is greater for the solid body than for the patch, which is consistent with Zong and Nepf's [29] results. This magnitude difference becomes smaller moving further downstream of the patch.

**Figure 12.** Spectral densities of the transverse velocity *Sv* located at (**a**) C5 and (**b**) C11 for Case SVF5 and Solid Case.

Figure 12b shows the change in longitudinal velocity spectral density with frequency for Case SVF5 and the Solid Case. When the frequency exceeds 0.16 Hz for the Solid Case and 0.04 Hz for SVF5, the linear decrease of spectral densities is in agreement with the Kolmogorov line, with a slope of −5/3, as shown in Figure 12b.

Channel blockage D/B = 0.25 (Figure 3) in this study may have some effect on the vortex shedding frequency in our study. Chen et al. [35] and Coutanceau and Bouard [36] stated that both Rec and the Strouhal number (StD = *f*D/*Udm*) increase with increasing channel blockage. Here, Rec is the critical Reynolds number at which unbounded steady flow past a circular cylinder becomes unstable. Sahin and Owens [37] found that Rec and StD were 285 and 0.44, respectively, for a blockage of 0.64. Zong and Nepf [29] considered D/B = 0.10, 0.18, and 0.35 for the emergent patch in their study and showed that ReD (=*Udm*D/ν) is *O* (104), which is within the turbulent wake regime. In our study, ReD is *O* (105), which is the within the turbulent wake regime. As a result, the channel blockage might have some effects on the vortex shedding frequency in our study.

Rominger and Nepf [31] defined deceleration through a long porous patch where length is much larger than width. They mentioned that the velocity at the centerline starts to decline from some distance upstream of the patch, which was also observed in this study (Figure 7a), and that the deceleration continues into the patch over a length xD, after which the flow reaches a steady velocity, U0—the final interior velocity. Rominger and Nepf [31] obtained that the interior adjustment length (xD) and the final interior velocity (U0) depend on the dimensionless parameter, *CDaD*, which is called the patch flow-blockage. Here, CD is the drag coefficient for the cylinders within the patch.

A low blockage effect is present for *CDaD* < 4, where interior velocity is driven by turbulent stress; while in a high blockage effect, when *CDaD* > 4, interior velocity is driven by a pressure gradient [29,31]. In our study, *CDaD* = 1.8 for φ = 0.05 (SVF5), *CDaD* = 7.5 for φ = 0.20 (SVF20), and *CDaD* = 13.8 for φ = 0.36 (SVF36), where CD = 1 for simplicity, as Zong and Nepf [29] considered in their study. There is a low blockage regime for the case φ = 0.05 < 0.1 and at this regime xD ~ 2 (CDa)<sup>−</sup><sup>1</sup> [31] and

xD = 33.33 cm in this study. Here, CDa is the canopy drag length scale, which is the length scale of flow deceleration associated with canopy drag. There is high blockage regime for the cases φ = 0.20 > 0.1 and φ = 0.36 > 0.1, which is consistent with Zong and Nepf's [29] experimental conditions. In these cases, the interior adjustment length xD ~D = 30 cm and diameter of the canopy patch provide us this length. As such, velocity just downstream of the patch is expected to be equal to the interior velocity for these cases.

#### **4. Conclusions**

The purpose of this study was to answer how the submergence ratio and density of vegetation patches affects the flow structure along the axis of the patch, downstream of the suspended cylindrical vegetation patch. For this purpose, the velocity measurements with ADV (Acoustic Doppler Velocimeter) were performed along the flow depth at the axis of the patch along different longitudinal distances downstream of the suspended patch for different submergence ratios (hv1/H = −0.33 and hv2/H = 0.67) and different patch porosities (φ<sup>1</sup> = 0.05, SVF5; φ<sup>2</sup> = 0.20, SVF20; φ<sup>3</sup> = 0.36, SVF36; φ<sup>4</sup> = 1.00, Solid Case).

It was found that there is a recirculation zone behind the obstruction, which is strongest for the highest density (SVF36) patch, while there is no recirculation zone for the lowest density (SVF5). It was also shown that U velocity increases more rapidly to the *Udm* value at the highest density (SVF36) compared to the lower densities (SVF5 and SVF20). Moreover, a steady wake region was observed downstream of the patch for SVF5, which is not present at higher SVFs.

The flow becomes stronger downstream of the gap region with increasing submergence ratio and flow behavior becomes similar to wall jet flow. It was found that SVF has an important impact on this flow structure in terms of the magnitude and the placement of the maximum velocity.

The patch with the lower depth ratio (hv/H = 0.33) affected the flow slightly downstream for the solid body and lowest SVF. However, suspended vegetation still affects the flow downstream of the gap region for a higher density of patches (SVF20 and SVF36) and it is concluded that high densities of submerged vegetation patches still have an influence on the flow structure downstream of the gap, even at low submergence ratios.

It was found that suspended vegetation patches do not have a significant effect on transversal velocity downstream of the patch for a low solid volume fraction (SVF5) at both depth ratios. As the suspended canopy becomes denser, the magnitude of transversal velocity increases, and the placement of diverging flow occurs closer to the patch, in the downstream section, with increasing SVF.

The effect of suspended vegetation on turbulence starts in the upstream section and causes a peak upstream of the gap region. It was found that the shear layer caused by the patch increases the turbulence level beneath the patch. It was also found that the individual stem turbulence is not significant for the highest patch density (SVF36), whereas it is significant at lower patch densities (SVF5 and SVF20). Individual cylinders cause the small-scale turbulence, which results in a higher turbulence region behind the patch. Moreover, the turbulence level decreases with increasing patch density (SVF36) and this turbulence level is close to the turbulence level for the solid body at mid-depth. This means that the suspended patch with the highest porosity patch (SVF36) behaves like a solid body in terms of turbulence and structure of the flow.

**Author Contributions:** A.Y.O. and D.Y. conceived and designed the experiments; A.Y.O. and D.Y. performed the experiments; A.Y.O. and D.Y. analyzed the data; A.Y.O. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Aydın Adnan Menderes University, Civil Engineering Department students Omer Botan YILDIZHAN, Ali CAMBOLAT and Ebru ERSIN who helped during the experiments.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
