A Two‐Layer Turbulence‐Based Model to Predict Suspended Sediment Concentration in Flows With Aquatic Vegetation

Traditional bed shear stress‐based models (e.g., Rouse model) derived from the classic parabolic profile of eddy viscosity in open‐channel flows fail to accurately predict suspended sediment concentration (SSC) in flows with aquatic vegetation. We developed a two‐layer, turbulence‐based model to predict SSC profiles in emergent vegetated flows. Turbulence generated from vegetation, bed, and coherent structures caused by stem‐bed‐flow interaction are considered into the near‐bed turbulent kinetic energy (TKE) to calculate the effective bed shear velocity, ubeff* . The model, validated by experimental data, further showed that the thickness height of the near‐bed layer (effective bottom boundary layer), Hb, varies with flow velocity and canopy density. Two additional models are provided to estimate Hb and ubeff* . The model is expected to provide critical information to future studies on sediment transport, landscape evolution, and water quality management in vegetated streams, wetlands, and estuaries.

transport and resuspension models relied on measurements of bed shear stress (e.g., Engelund & Fredsoe, 1976;Meyer-Peter & Muller, 1948;Nino & Garcia, 1998;Parker, 1979;van Rijn, 1984;Yalin, 1963). However, studies have shown that bed shear stress models do not work in regions with vegetation (Tinoco & Coco, 2014;Yager & Schmeeckle, 2013). Those models do not account for the turbulence generated by vegetation (Nepf, 2012), which has been recognized as a prominent feature in natural flow environments that drives sediment resuspension Sumer et al., 2003). Through laboratory experiments, Tinoco and Coco (2016) have demonstrated how near-bed vegetation-generated turbulence affects the incipient motion of sediment in vegetated flows. Yang and Nepf (2018) proposed a turbulence-based bed-load transport model based on the relation between turbulence-induced lift force and turbulent kinetic energy (TKE), similar to the nonvegetated case discussed by Nino and Garcia (1996). Their model showed a good agreement with experimental data when including both TKE components from vegetative turbulence and bed-shear turbulence.
To construct a more accurate model for predicting SSC, we must consider the local flux of sediment entrained from the bed by near-bed coherent structures, generated by stem-bed-flow interaction (e.g., horseshoe vortex). Previous studies have revealed that these coherent structures can significantly alter the near-bed stress locally in front of cylinders (Schanderl et al., 2017a(Schanderl et al., , 2017b, resulting in large differences between the time-averaged and instantaneous turbulent quantities (Apsilidis et al., 2015). These fluctuating coherent structures can cause nonnegligible contributions to sediment resuspension (Cheng et al., 2018;Li et al., 2018). In response to the above effects, we propose a new framework considering all the turbulence components from (1) vegetation, (2) bed shear, and (3) fluctuating coherent structures, into nearbed TKE calculation for modeling SSC in flows through emergent vegetation. The new turbulence-based SSC model is based on a two-layer effective turbulent stress profile, and accurately predicts SSC measured experimentally.

Two-Layer SSC Model
Starting from a traditional approach that considers sediment particles as continuous diffusive scalars, based on equilibrium of suspended sediment without reacting with the fluid, turbulence stress is balanced with sediment's settling as: where C is the time-averaged SSC, w s is the sediment settling velocity, and D t is the flow eddy diffusivity in the vertical direction. D t is usually equated with the vertical flow eddy viscosity, ν t , assuming the turbulent Schmidt number, Sc t = ν t / D t = 1, according to Reynolds analogy of eddy mixing. By solving the above differential equation, we can obtain the general form of the SSC profile where C b is the time-averaged SSC at a reference z-location, b, near the bed. Rouse profile of SSC (Rouse, 1937) has been widely applied in open-channel flows based on a parabolic profile of eddy viscosity, ν t , as: u is the bed shear velocity and κ is the von Karman constant = 0.41. Replacing Equation 3 into Equation 2, the Rouse profile of SSC can be obtained as: However, in vegetated flows, the profile of ν t is no longer parabolic (Defina & Bixio, 2005), thus Rouse profile cannot be used to accurately predict SSC (Lopez & Garcia, 1998).
According to the turbulent-viscosity hypothesis, ν t can be obtained by the ratio between the turbulent stress, τ t , and the mean flow velocity gradient, dU/dz, which has been proved valid for simple shear flows (e.g., open channel flows). However, in canopy flows, counter-gradient momentum transport is generally observed close to the bed (Shaw, 1977) due to stem-bed-flow interaction, which could modify the turbulent-viscosity hypothesis in the near-bed region (Raupach & Thom, 1981). Hence, we introduce a new, more effective approximation for this near-bed region. Above the near-bed region, Equation 6 still holds and suggests that profiles of mean flow velocity and turbulent stress are key to derive ν t in vegetated channels. Previous studies have shown a characteristic constant velocity profile above a certain distance from the bed in emergent vegetated flows (Nepf & Vivoni, 2000;D. Liu et al., 2008;Stoesser et al., 2010). Yang et al. (2015) proposed a two-layer linear stress model to predict bed shear stress and the velocity profile in vegetated flows under smooth bed environments. Etminan et al. (2018) further provided a balance argument between TKE production and TKE dissipation rate to precisely estimate the thickness of the viscous layer. However, most of the bottom boundaries in natural water environments are hydraulically rough, where roughness elements (e.g., sediment particles) penetrate through the viscous layer. When k s /δ ν > 1, where k s is the bottom roughness height, and δ ν is the height of viscous layer, viscous layer becomes relatively small and thus negligible. Also, varying bottom roughness can further increase the boundary layer thickness (Hanmaiahgari et al., 2017;Schultz & Flack, 2007).
To account for the bottom roughness effect, we propose a modified two-layer SSC model based on an effective turbulent stress profile for flows with emergent vegetation. In the region close to the bed but above the viscous layer, turbulence generated from both the bed and the vegetation is dominant, which causes counter-gradient momentum transport due to interacting stem-bed-flow coherent structures. To account for such complex interactions, we approximated the turbulent stress profile as an effective linear-stress profile, following Yang et al. (2015), and included the local flux of sediment entrained from the bed due to fluctuating near-bed coherent structures into the near-bed TKE calculation to update the linear turbulent stress (see details in Section 4.4). This region is so called the linear-stress layer or the effective bottom boundary layer. Above the linear-stress layer, the velocity profile reaches a constant value, and the turbulent stress in this region also approaches a constant (      . t u w const ). This region is called the constant velocity layer or the constant stress layer. In this layer, turbulence generated from the bed can be ignored, and sediment suspension is mainly supported by the stem-generated turbulence.
By setting H b as the thickness height of the effective bottom boundary layer, the effective turbulent stress can be expressed as: where H is the water depth, and * b eff u is the effective bed shear velocity, derived from the total near-bed turbulent stress, which considers turbulence generated from the bed shear stress, vegetation, and the coherent structures. * b eff u is different from the bed shear velocity, * b u , in Equation 3, which was derived directly from the bed-shear stress or bed-shear turbulence. More details are described in Section 4.4. Since log-law is assumed valid for the effective velocity profile in the linear stress layer, the effective mean flow velocity gradient can be estimated as Replacing Equations 7 and 8 into Equation 6, the reciprocal of the effective eddy viscosity, ν eff , of the two-layer model can be obtained Finally, by further replacing Equation 9 into Equation 2 and solving the integral of 1/ν eff , we can obtain the final form of the two-layer model for SSC with an updated Rouse number, In the above model, there are two critical parameters that need to be determined: the thickness height of the bottom boundary layer, H b , and the effective bed shear velocity,  b eff u . We developed two additional models to estimate these parameters in Sections 4.3 and 4.4, respectively.

Experimental Setup
The experiment was conducted in an Odell-Kovasznay type recirculating flume (Odell & Kovasznay, 1971, Figure 1a), with a straight test section 2 m long, 0.15 m wide, and 0.6 m deep. A vertical axis disk-pump with uniformly distributed disks drives the flow to produce a uniform velocity profile with minimal vertical disturbance. As water continues to recirculate in an essentially infinite loop, the design guarantees full development of the boundary layer and turbulence features in the flume, which ensures applicability of Equation 1 with an assumption that the sediment profile reached equilibrium and fully developed with the flow. The disk rotational frequency, ω, is set by the frequency, f, of an inverter as ω (rpm) = 6.6f (Hz).
High-stiffness aquatic vegetation such as Sagittaria sagittifolia and Sparganium erectum are mimicked by an array of rigid acrylic cylinders with diameter d = 0.64 cm. The water depth, H, is set as 10 cm. The length of vegetation array is 156.2 cm (246 d), and the leading edge is 20.3 cm (32 d) downstream from the beginning of the straight test section. A 7.6 cm (12 d) gap was left within the array from 135.3 cm to 142.9 cm (213 d-225 d) from the beginning of the array for Particle Image Velocimetry (PIV) and SSC measurement. A staggered configuration was used with two array densities to cover from sparse to dense conditions, ah = (0.1-0.5) (Nepf, 2012). The detailed array configuration is shown in Figure 1b.
To account for scaling of sediment dynamics in laboratory experiments (Heller, 2011;Kleinhans et al., 2014), crushed walnut shells are chosen as sediment substrate (Briskin et al., 2002), with a 10 cm thick sediment bed as the initial setup in the straight test region. The sediment density is ρ s = 1.2 g/cm 3 (water density, ρ = 1.0 g/cm 3 ) with a diameter range 0.84-1.68 mm (median grain size, D s ≈ 1 mm). The walnut shell settling velocity is w s ≈ 1 cm/s, which is approximated by Rubey (1933): To ensure sediment resuspension in our tests, mean flow velocities ranged from 1.7 cm/s to 6 cm/s within the vegetation canopy (f = 20-40 Hz). This range exceeds the estimated critical velocities (U c ≈ 1.7 cm/s) and 1.5 cm/s under sparse (ah = 0.1) and dense (ah = 0.5) array (with the light-weight sediment, respectively). Critical velocities were estimated by first obtaining the critical velocity under bare-bed condition by Ikeda (1982): where C L = 0.85C D , F = f (U 0 ), and then updating the values for the vegetated flows based on Yang et al. (2016): where U c is the critical velocity in vegetated flows, U 0 is the critical velocity under bare-bed condition, ϕ is the solid volume fraction of the canopy (≈(π/4)ad for circular cylinders), B is the coefficient noted as C in Yang et al. (2016).
Five experimental runs with emergent vegetation arrays and three test runs with bare-bed setup were conducted (Table 1). Under the same f, the mean flow velocity was much higher in the bare-bed cases because of the extra drag exerted by the vegetation. However, the turbulence generated by the vegetation enhanced resuspension, which made the number of suspended particles with vegetation at low speeds similar to the bare-bed TSENG AND TINOCO 10.1029/2020GL091255 5 of 14  case at high speeds. Maximum velocities were capped to ensure adequate capture of individual suspended sediment particles with our quantitative imaging system (details in Section 3.3).

Hydrodynamics Measurements
A 5 W continuous-wave laser system was used to generate a vertical planar light sheet for PIV measurement with a thickness of < 1 mm at the centerline of the flume. The same setup was used to generate a horizontal light sheet 2 cm below the surface to investigate the horizontal flow structure and validate the 2D vertical-plane approach of the study in our thin flume, where wall effects can be ignored, even for the sparse and bare-bed cases (see details in the supporting information). The light sheet covers the whole observation gap (12 d long) (red-dashed box in Figure 1 (c)). This gap was left to avoid laser reflection issues by the acrylic cylinders within the array. As shown in the supporting information ( Figure S1), velocity measured within the gap is representative of velocity within the array. A 5-Megapixel CCD camera, JAI GO-5000M-USB3, with a Navitar 25 mm focal length was used to capture 8 -bit grayscale images at 60 Hz for 1 min (3,600 images) in each run. Raw images were processed in PIVlab (Thielicke & Stamhuis, 2014). The image spatial calibration factor, CF, for each case is shown in Table 1. Three consecutive 50% size passes with 50% overlapped interrogation areas were used to obtain higher resolution results during cross correlation calculations. The sub-window size of the first pass followed the one-quarter rule (Keane & Adrian, 1990) (detailed size information in Table 1).

SSC Measurements
SSC was estimated by counting the time-averaged particle pixels in the red-dashed observation area in Figure 1(c), to quantify the probability of sediment particles to be present at each location. We used bright lighting with white background to contrast particles (dark pixels) and ambient water (white pixels, Figure 1(d)). Using the same camera as for PIV, 150 images were taken at 0.5 Hz frame rate for each case. We set a threshold intensity to identify sediment particles. We specify a binary value to each pixel (1 = sediment, 0 = water), and average over all frames and along the x-direction to get a time-averaged SSC profile, C(z), as a probability of the presence of sediment grains in the water column. The bed location, z = 0, is determined by the top location where C(z) = 100%. Calibration for SSC was done by adding known amount of sediments into a small tank with known volume of water, continuously stirred to ensure that sediment particles are well-mixed and uniformly distributed. Following the same method (intensity threshold, binary images, temporal and spatial averaging) with a 5-minute-long record for each known concentration. The near-bed concentration, C b , for each case is provided in Table 1.
The quantitative imaging method for SSC works well in our thin flume. Sediment grain size was larger than the image pixel size to avoid overestimation of the measured SSC. The amount of suspended sediment particles was limited to an appropriate range to avoid saturation, which would cause extremely dense clouds of suspended sediment that would prevent us from getting an accurate metric. In low-speed cases, when there are not many suspended particles, the camera resolution and the size of sediment particles determine a lower-limit baseline, which occurs at some specific z-locations where just a few grains appear on a small number of frames. This is more noticeable in the upper water column near the free surface.
To improve the analysis at low SSC (i.e., low number of recorded particles), a moving-average-frame technique was applied to include the adjacent concentration data into the analysis. The size of the moving frame was designed as 3-5 times the sediment size to ensure that the averaging SSC data at each z-location included the entire sediment particle, and averaged over the adjacent region within twice the particle size. More details are described in the supporting information.

SSC Profiles of Bare-Bed Channels
To validate the SSC observation method, three bare-bed test cases under different mean flow velocities were conducted and compared with the classical Rouse profile (Equation 4) as shown in Figure 2. The near-bed TSENG AND TINOCO 10.1029/2020GL091255 reference location for the Rouse profile was set as b = D s ≈ 1 mm. When the flow velocity goes higher, more sediments are suspended, which makes higher SSC in the water column, showing better agreement between the experimental data and the fitted Rouse profile as R 2 value is as high as 0.86 (Figure 2(c)). Under lower flow velocity, there are fewer sediments suspended in the water column, resulting in lower SSC profile and larger deviations from the fitted Rouse profile as R 2 value drops to 0.65 (Figure 2(a)). The deviation, due to the camera resolution and the particle size, is discussed in Section 3.3. However, the fitting coefficient still shows a satisfactory fit, ensuring R 2 > 0.65 among all the cases. In Figure 2

SSC Profiles of Emergent Vegetated Channels
By applying the proposed two-layer model derived in Section 2, Figure 3 shows the model fitted results against the measured SSC data points within emergent arrays. Like the bare-bed case, the near-bed reference location for the model is also set as b = D s ≈ 1 mm. The height of the boundary layer region is defined by the vertical location where mean velocity profile reaches a constant value as shown by the dotted line in Figure 3. Within the boundary layer region, flow structure is strongly affected by the flow interaction between bed shear and vegetation stems. The near-bed velocity bulge is one of the features caused by flow coherent structures (e.g., horseshoe vortex and junction vortex), which entrains fast moving fluid with high momentum from the surrounding region into the wake (D. Liu et al., 2008;Stoesser et al., 2010). As a consequence, the varying velocity profile in this region alters SSC, while above H b , the unchanged mean flow velocity makes a constant SSC profile. The model performs well to predict SSC throughout the entire water TSENG AND TINOCO 10.1029/2020GL091255 8 of 14 with mean velocity and array density, such that their parameterization is needed for accurate predictions.

A Model to Predict Thickness Height of the Effective Bottom Boundary Layer
According to previous studies, within the effective bottom boundary layer, the layer thickness height, H b , is governed by an equilibrium between rates of turbulence production, P, and the viscous dissipation, ϵ, (Brouwers, 2007;Karimpour & Venayagamoorthy, 2013): Following a similar argument proposed by Etminan et al. (2018) that TKE production in the cylinder wakes balances the rate of work done by form drag (Tanino & Nepf, 2008), turbulence production, P, can thus be expressed as the wake production, P w , generated by the vegetation. Equation 14 then becomes P w can be further scaled as (Nepf, 1999;Nepf & Vivoni, 2000): where C D is the drag coefficient of the vegetation canopy, and ϵ can be scaled as: Replacing Equations 16 and 17 back into Equation 15, we obtain: In our case, since the roughness density of the canopy is in a range of ah = 0.1 to 0.5, C D ≈ 1 can be a valid approximation (Nepf, 1999). Etminan et al. (2018) followed the linear stress model by Yang et al. (2015) to propose an approximation of the bed shear velocity.
where k s is the bottom roughness. The above assumption of the inverse proportion between  b u and k s is only valid when k s /H ≪ 1. k s is usually estimated from the sediment size (k s = D s ≈ 1 mm in our cases). Finally, replacing Equation 20 into Equation 18, we estimate H b based on bottom roughness as: where α is an empirical coefficient that needs to be experimentally determined. The new proposed H b estimation in Equation 21 accounts for the bottom roughness effect in natural water environments.
The model suggests that H b increases as U or a decrease, meaning that when the flow velocity is low, the effective bottom boundary layer is thicker, similarly to the turbulent boundary layer becoming thicker when Reynolds number is lower. Also, for sparser vegetation, the flow can behave more like an open-channel flow, with a larger portion of the flow following the logarithmic profile. Higher k s will also increase H b because stronger turbulence can be generated from the rougher bed. However, the model is asymptotic to infinity when a 2 UH 4 /νk s → 0 as shown in Figure 4(b), requiring a lower limit of this model according to H b /H ≤ 1, which yields: Yang and Nepf (2018) showed that sediment incipient motion in vegetated channels depends on the total near-bed TKE, k t , which is assumed to be the sum of bed-generated turbulence, k tb , and vegetation-generated turbulence, k tv . Then k tb can be approximated by (Stapleton & Huntley, 1995):

A Turbulence-Based Model to Predict Effective Bed Shear Velocity
where C ub is the integral coefficient that converts mean flow velocity to bed shear velocity based on log-law (Wilcock, 1996):  value, c, varies under different Re d and ah. High SSC values can potentially cause turbulence modulation in the flow. However, the ranges of SSC observed in our study are low enough to avoid such particle feedbacks and potential turbulence stratification effects (Cao et al., 2003). High SSC and a sharp SSC gradient only occur in the near-bed region, which could introduce particle feedbacks to influence the magnitude of c. To obtain a more comprehensive range of c would require more experimental or numerical studies, considering a broad range of sediment size, density, and vegetation characteristics, which is out of the scope of the current study.

Model Application
According to the proposed two-layer model and the associated turbulence-based parameterizations of H b and  b eff u , we can predict SSC in flows with emergent vegetation by following the procedure summarized as follows:

Conclusions
The vegetation-generated turbulence alters the flow structure and interacts with bed-generated turbulence, creating near-bed coherent structures that enhance sediment suspension. The traditional Rouse profile derived from the classic parabolic profile of ν t in open-channel flows fails to accurately predict SSC in regions with aquatic vegetation. We propose a two-layer model based on mean flow and turbulent stress profiles with two additional turbulence-based models that provide estimations of the effective bottom boundary layer thickness, H b , and the effective bed shear velocity  b eff u . The models show that H b increases as U or a decrease, while  b eff u is determined by a combination of k tv and k tb with an empirical constant, c, to parameterize the local fluxes of sediment entrainment by near-bed coherent structures. More experimental or numerical studies are needed to obtain a more comprehensive expression of c, which varies under different Re d and ah. Experimental results and observations from previous studies validated the newly developed two-layer model under ah = {0.1-0.5} and Re d = {100-400}, which provides better prediction and more detailed insights on sediment resuspension and SSC in emergent vegetated flows. The model is expected to deliver critical information for future studies on landscape evolution and water quality management in vegetated channels, wetlands, and estuaries.