Effect of submerged vegetation on solute transport in an open channel using large eddy simulation

Advances in Water Resources 97 (2016) 87–99

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Effect of submerged vegetation on solute transport in an open channel using large eddy simulation J Lu∗, HC Dai College of Water Conservancy and Hydropower, Hohai University, Xikang Road., Nanjing 210098, PR China

a r t i c l e

i n f o

Article history: Received 17 March 2016 Revised 20 August 2016 Accepted 2 September 2016 Available online 5 September 2016 Keywords: LES Vegetation Scalar transport Quadrant analysis

a b s t r a c t Existence of vegetation plays a significant effect on the flow velocity distributions, turbulence structures and solute mixing in an open channel. This paper has implemented a 3D large eddy simulation model for the flow and scalar transport in the open channel with vegetation. The model can produce a typical turbulence characteristics and concentration distribution with vegetation. The scalar transport mechanism is quantitatively explained by the turbulent Schmidt number, Reynolds flux, coherent structures and quadrant conditional analysis. A dominance of ejection-sweeping events occurs in the process of the momentum and scalar flux transport. The spectral analysis is used to identify the Kelvin–Helmholtz frequency. The turbulence characteristics of the length scale of vortexes, Kelvin–Helmholtz frequency and Reynolds stress etc. are analyzed with the vegetation density. The model quantitatively predicts the trend of decreasing in the concentration distribution along the flow direction with the increasing of vegetation density. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Vegetation plays a significant role in the sustainable development in streams and rivers, and it has important effects on the flow characteristics such as the velocity distributions, turbulence structures and the process of solute mixing (Tsujimoto, 1999; Nepf and Ghisalberti, 2008; Okamoto and Nezu, 2010). As a result, it is important for environmental or ecological engineers to investigate the hydrodynamic characteristics and process of solute transport through submerged vegetation. Effects of vegetation on the flow characteristics have been studied by a lot of laboratory flume experiment works (e.g. Carollo et al., 2002; Järvelä, 2002; Huthoff et al., 2007; Poggi et al., 2004). The mean velocity distribution within submerged vegetation does not follow the universal logarithmic law. Moreover, the peak value of the Reynolds stress occurs at the canopy edge and decreases linearly above the canopy to a value of zero at the free surface. In a study involving the submerged vegetation study, Nezu and Sanjou (2008) found that the entire vertical velocity distribution can be divided into three layers: emergent zone, mixing-layer zone and log-law zone (e.g. Neary, 2003; Huai et al., 2009). The authors divided the zone from the river bed to the 90% of the maximum value of the stress as the emergent zone, in which turbulence is dominated by stem-scale turbulence. They defined the whole range from the elevation of the mean velocity starting to obey the ∗

Corresponding author. E-mail address: [email protected] (J. Lu).

http://dx.doi.org/10.1016/j.advwatres.2016.09.003 0309-1708/© 2016 Elsevier Ltd. All rights reserved.

log-law to the water surface as the log-law zone. The region between the emergent and log-law zone is the mixing-layer zone, in which turbulence is governed by the KH vortices. The difference in drag magnitude between the non-vegetated and vegetated zones leads to Kelvin–Helmholtz (KH) vortices occurring at the interface between the non-vegetated and vegetated zones. The KH vortices can promote the mass and momentum transport between the within and over canopies (Raupach et al., 1996; Nepf, 1999; Nepf and Vivoni, 20 0 0). The KH instabilities have significant effects on the large-scale turbulence structures and the momentum transfer between the non-vegetated and vegetated regimes. The large-scale turbulence structure dominates the turbulent diffusion in the mixing-layer zone and the turbulent diffusion coefficient is not unity (Ghisalberti and Nepf, 2005). The above works provide physical insight of the flow phenomenon with vegetation. However, many numerical models have been recently developed for free surface flows with vegetation. For examples, Li and Yan (2007) and Busari and Li (2015) employed the resistance model with the Spalart–Allmaras models (Spalart et al., 1997), in which the presence of vegetation is obtained in the turbulence model by introducing the drag-related source terms, to simulate the interactions between the flows and vegetation. The two-equation resistance models based on the Reynolds-averaged Navier–Stokes (RANS) models was proposed and validated (e.g. Shimizu and Tsujimoto, 1994; Lopez and Garcia, 1998; Leu et al., 2008). The anisotropic Reynolds stress models (RSM) for flows with vegetation were developed by Naot et al. (1996) and Choi and Kang (2004). Moreover, a large eddy simulation, which not only predicts

88

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

the time-averaged velocity fields but also can give instantaneous velocity fields and resolved turbulence structures, has been used for analyzing the instantaneous monument transferring process at the interface between the vegetated and non-vegetated region (Su and Li, 2002; Cui and Neary, 2008; Li and Xie, 2011; Lu and Dai, 2016a). In addition, the presence of vegetation can alter the flow of streams and rivers by affecting the bed shear stress, sediment concentration and water quality and so on (Schultz, 2002; Nepf et al., 1997). The vegetation can also affect the solute diffusion due to the blocking of vegetation resulting in the mechanical dispersion phenomenon (Nepf, 2012). The mathematical models considering the mechanical dispersion were recently developed by Zhang et al. (2010), Poggi et al. (2006) and Lu and Dai (2016b). Moreover, a series of scalar flux models were proposed by Kang and Choi (2009). The authors investigated the averaged concentration profiles deviating from the Gaussian distribution in a submerged vegetation study. Okamoto and Nezu (2010) employed a LES model to find that the vertical turbulent exchange contributes largely to the scalar transport in the mixing layer zone. Although extensive research works mentioned above on the flow characteristics in a vegetated channel and a few studies on the solute mixing process in an open channel with vegetation have been investigated, the process of the solute transport within the vegetation has not been fully investigated. In the present work, a large eddy simulation model was used to study the effects of vegetation on the solute transport in an open vegetated channel. Therefore, the flow and solute concentration characteristics were investigated. The effects of submerged vegetation on the concentration transport were analyzed. The quadrant and spectral analyses were used to understand the flow structure and solute transport mechanics. 2. Governing equations The basic equations consist of the filtered continuity equation, N-S equations including additional drag terms produced by vegetation which is evaluated as an internal source of resistant drag force per unit fluid mass, and the concentration transport equation. The governing equations for incompressible fluid flow in the Cartesian coordinate can be written in following (Lu and Dai, 2016b):

∂ ( u ) = 0, ∂ xi i

(1)

  1 ∂p 1 ∂ ( ui ) ∂ ∂ ∂u + ( u i u j ) = gi − + υ i − τi j − Fi , (2) ∂t ∂xj ρ ∂ xi ∂ x j ∂xj ρ   ∂C ∂ ∂ ∂C + (u C ) = D + Sc , ∂t ∂ x j j ∂x j t ∂xj

(3)

where ui are the filtered velocity components in the Cartesian coordinate, xi are the Cartesian coordinate, t is the time, gi are the gravity force components, p is the pressure, ρ is the density of fluid, υ is the molecular viscosity, Fi are the resistance force components per unit volume induced by vegetation in the x, y and z directions, respectively. C is the concentration, Dt is the diffusion coefficient and Sc is the source term . The sub-grid stresses τi j can be written as:

1 3

τi j − τkk δi j = −2υt Si j , Si j =

1 2



∂ ui ∂ u j + ∂ x j ∂ xi



(4)

Here, the eddy viscosity within vegetation was modeled by the maximum value of the sub-grid scale viscosity (υsmg ) and the eddy viscosity generated by the stems (υveg ), which was proposed by Lu

and Dai (2016a).

υt = max (υveg , υsmg )   = max

 1 2 Cu Cd DN (ui ) u j u j 2

1 / 2

 min (D, S ), (Cs ) |S| , (5) 2

where Cd is the drag coefficient of stem, D is the vegetation stem diameter or width normal to the flow direction, N is the number density and Cs is the Smargorinky constant.  = (xyz )1/3 , 1/2 |S| = (2Si j Si j ) , Cu =0.09. The diffusion coefficient Dt within vegetation considering the mechanical dispersion (Lu and Dai, 2016b) is given as:



Dt =

v+

β2 2



N um D3 + υt /St ,

(6)

√ where β = 2, St =1 is the turbulent Schmidt number and um is the absolute velocity within vegetation. The first and second terms of right hand in the Eq. (6) are the modular and mechanical diffusion terms, respectively. The third term of right hand in the Eq. (6) is the turbulent diffusion, which could be replaced by the Eq. (5). The average force per unit volume within the vegetation domain is given as follows:

Fi =

  1 1 Cd ρ DN ui u j u j = Cd ρα ui u j u j , 2 2

(7)

where α is the vegetation density. 3. Numerical methods and boundary conditions Since the free surface elevation varies with time, and the bottom is uneven for hydrodynamic flows, it is difficult to discretize the domain along the vertical direction. To solve the uneven physical domain, the governing equations are first transformed into the vertical σ -coordinate. The splitting operator approach is then implemented to numerically solve the governing equations (Lu and Dai, 2016a). At each time interval the N-S equations are solved in three steps: advection, diffusion (source/sink) and pressure propagation steps, and the concentration transport equation is solved in two steps: advection and diffusion steps. The combination of Lax– Wendroff method and quadratic backward characteristic method on non-uniform grids is used to solve the advection equations. The central difference scheme in space is used to solve the diffusion and pressure propagation equations. A bi-conjugate gradient stabilized method (BI-CGSTAB) is used in this paper to solve the Poisson equation (Lu et al., 2010). For solving the governing equations, boundary conditions must be specified. Different boundary conditions are usually used in free surface flow problems. At the free surface the dynamic and kinematic condition is applied, which is written as:

∂η ∂η ∂η = u3 − u1 − u2 , ∂t ∂x ∂y

(8)

in which η is the surface elevation. At the inflow boundary, the convective boundary condition for the surface elevation is imposed, and the mean velocity components are specified, and the fluctuating velocity proposed by Jarrin et al. (2006) is used.

∂η ∂η = −Uc , u = Ui + ui , ∂t ∂x i 

(9)

where Uc = gh, Ui is the mean velocity components. The fluctuating velocity components ui is generated by the synthetic eddy method of Jarrin et al. (2006). The synthetic eddy method based on the characterization of turbulence as a series of coherent eddies which has been implemented in our previous research (Lu and Dai, 2016a), is used here.

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

89

Table 1 Computed parameter of cases. Run Run Run Run Run Run Run Run Run Run Run Run Run Run Run Run

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Vegetation height(m)

Vegetation density(m−1 )

Parameter (Cd α hv)

Solute location(y(m); z(m)

Water depth(m)

Bulk velocity(m/s)

/ 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05

/ 7.625 7.625 7.625 2 5 10 15 20 30 40 50 60 70 80

/ 0.38 0.38 0.38 0.1 0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4

0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2; 0.2;

0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15

0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12

0.05 0.05 0.025 0.075 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05

Note: We used Cd = 1.0 suggest by Schultz (2002) (Nepf (2012)) to calculate the parameter of Cd α hv.

The convective boundary condition is imposed at the outflow boundary.

∂φ ∂φ = −Uc , ∂t ∂x

(10)



where φ = η or ui . When φ = η, Uc = gh; when φ = ui , Uc is the convective velocity at the outflow boundary. The no-slip boundary condition for velocity is used on the solid walls. The wall function proposed by Werner and Wengle (1993) is used on the solid wall.

|τub | =

2μ|u p | n ,

|τub | = ρ +



1+B A





|u p | ≤

μ

2 ρ n

A2/(1−B)

μ 1+B 1−B [(1+B )/(1−B )] A 2 ρ n



× ρμ n

B

2/(1+B) , |u p |

,

|u p | >

μ

2 ρ n

(11)

A2/(1−B)

where A = 8.3, B = 1/7, τ ub is the wall shear stress, up is the velocity of the first interior point, n is the distance to the solid wall and μ is the dynamic viscosity coefficient. For concentration equation, the zero gradient is set to the normal direction at the boundaries for simplicity, namely,

∂C = 0, ∂n

(12)

in which n is the normal direction.

(2010), was used to validate this study. In their experiments, the rectangular flume was of length 10.0 m, width 0.40 m, and depth 0.30 m. Vegetation was simulated by rigid-strip plates with 8 mm width, 1 mm thickness and 50 mm length. A continuous source of tracer with initial concentration C0 was released at tip of vegetation and the injection velocity was kept at the local velocity. A particle image velocimetry (PIV) and laser-induced fluorescence (LIF) was used to measure the flow and concentration fields. The case used for comparisons has the following parameters: the flow depth (H) is 0.15 m, the bulk velocity (U) is 0.12 m/s and the vegetation density (α ) is 7.625 m−1 (Referred to Run 2 in Table 1). The drag coefficient of vegetation Cd = 1.10 is used. The computational domain is 15 m long, 0.4 m wide and 0.15 height. The non-uniform grids of about 0.65millilion are adopted, with the minimum grid size of 0.00375 m and the maximum grid size of 0.02 m. x, y and z are the streamwise, spanwise and vertical coordinates, respectively. The time step is 0.002 s. The time step is adaptive during the simulation, following the Courant–Friedrichs– Lewy (CFL) criterion, which is taken as 0.65 to ensure numerical stability.

5. Results and discussions 5.1. Mean velocity distribution

4. Problem setup The laboratory flume experiment on the solute transport in a vegetated channel, which was carried out by Okamoto and Nezu

Fig. 1. Computed vertical mean velocity profile for Run 2 and measured data by Okamoto and Nezu (2010).

Before investigation the concentration field, it is of importance to validate the mean velocity in the open channel with vegetation. Fig. 1 compares the computed normalized mean streamwise velocity at the mid-plane in the transverse direction for Run 2 and the experimental data by Okamoto and Nezu (2010). The values are normalized by the time-averaged velocity at the vegetation edge uh . It can be seen that the computed profile is generally in close agreement with the experimental measurements. Because of the blocking of vegetation, a velocity defect is produced within the vegetation region, as shown in Fig. 2(b)–(m), leading to the velocity profiles departing from the typical logarithmic distribution in the open channel without vegetation (Fig. 2(a)). Fig. 3 compares the velocity profiles with the different vegetation density in the fully development region. Obviously, the presence of vegetation has a great effect on the velocity distribution within the vegetation zone. One can observe that the increasing the vegetation density decreases the velocity within the canopy and correspondingly increases the velocity above the canopy. The presence of vegetation results in the velocity profile deviating from the typical logarithmic distribution. However, the velocity profile decreases slightly within the vegetation and increases slowly above the vegetation with the vegetation density α ≥ 30 m−1 in our

90

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

Fig. 2. Computed normalized velocity contour: (a) for Run 1, (b) for Run 2, (c) for Run 5, (d) for Run 6, (e) for Run 7, (f) for Run 8, (g) for Run 9, (h) for Run 10, (i) for Run 11, (j) for Run 12, (k) for Run 13, (l) for Run 14 and (m) for Run 15(Dashed zone is the vegetation zone).

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

91

Fig. 3. Computed normalized velocity profiles at the mid-plane in the transverse direction.

Fig. 5. Computed normalized scalar concentration contours (a) for Run 1, (b) for Run 3, (c) for Run 4 (Dashed line is the vegetation edge).

Fig. 4. Computed normalized scalar concentration contours at the section of (a) y/hv = 4 and (b) z/hv = 1 for Run 2 (Dashed line is the vegetation edge).

cases. In other words, the velocity is not highly sensitive to the vegetation density α ≥ 30 m−1 .

Fig. 6. Computed normalized scalar concentration profiles at x/hv = 1. 6 from the source.

5.2. Mean concentration distribution Fig. 4 shows the computed passive scalar concentration contours normalized by C0 at the sections of y/hv = 4 and z/hv = 1 for Run 2. It can be observed that the spanwise shape of concentration distribution keeps the Gaussian type. However, the vertical shape of concentration distribution deviates from the Gaussian shape after a distance from the source, duo to the large variation of the streamwise velocity near the tip of vegetation leading in a steep velocity gradient (Fig. 3) and the large difference in the concentration diffusivity inside and above the vegetation. The details making the difference in the concentration diffusivity will be discussed in the following part. By comparison of the spanwise and the vertical shape of scalar concentration distribution (Fig. 4(a) and (b)) and comparison of the vertical shapes of scalar concentration distribution with vegetation or without (Figs. 4(a) and 5(a)), one

can easily conclude that both the variation of streamwise velocity gradient and the difference in the concentration diffusivity are major factors leading in the vertical shape of scalar concentration distribution deviating from the Gaussian type, which can be quantitatively calibrated by the scalar concentration profiles (Fig. 6). Fig. 4(a) also reveals that the local maximum concentration values are located roughly along the vegetation edge, which is consistent with the experimental observation of Okamoto and Nezu (2010). Nezu and Sanjou (2008) divided the open channel flow with submerged vegetation into the three sub-zones: emergent zone, mixing-layer zone and log-law zone. For Run 2, we calculated that the range of the emergent zone is about 0∼0.035 m; the scope of the mixing-layer zone is about 0.035∼0.064 m; and the log-law zone is about 0.064∼0.15 m. Thus, the injected dye of Run 2 is located in the mixing-layer zone; the injected dye of Run 3 is settled

92

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

Fig. 7. Comparisons of the concentration decaying along the vegetation edge with the different vegetation density.

Fig. 9. Computed distributions of the turbulence intensity for Run 2: (a) urms /U, (b) vrms /U and (c) wrms /U(Dashed line is vegetation edge).

Fig. 8. Computed concentration profiles for Run 2 and measured data by Okamoto and Nezu (a) located at x/hv = 1.28 from the source and (b) located at x/hv = 1.96 from the source.

in the emergent zone; and the injected dye of Run 4 is situated in the log-law zone. Figs. 4(a) and 5(b) and (c) show the contours of concentration distribution in the three sub-zones, respectively. By comparing Figs. 4(a) and 5, one can easily observe that the vertical mixing and diffusion rate of scalar concentration across the vegetated region is much faster than across the non-vegetation region. This is also shown by the scalar concentration profiles (Fig. 6). That maybe be explained by the reasons of both the vegetation resistance resulting in a steep vertical velocity gradient in the canopy and the difference in the concentration diffusivity inside and outside the vegetation mentioned above. The large difference in velocity magnitude at the layer between the vegetation and nonvegetation zone formats a shear layer, leading to the enhancement of the turbulent scalar diffusion. Fig. 7 compares the influence of the vegetation density on the concentration decaying along the vegetation edge. It is apparent that the existence of vegetation has a crucial effect on the concentration decaying along the vegetation top, and the value of scalar concentration decreases with the increasing of the vegetation density. It is also of interest to note that the concentration value decreases slightly along the vegetation edge with the vegetation density α ≥ 30 m−1 . This is the reason of that the increasing in vegetation density α ≥ 30 m−1 will not remarkable decrease the velocity distribution within the vegetation and not striking increase the velocity gradient near the canopy edge discussed before. A turbulent Schmidt number is an empirical parameter of describing the pollution concentration diffusion and mixing. The turbulent Schmidt number is not unity in vegetated flows. The turbulent Schmidt number is rough 0.47 experimented by Ghisalberti and Nepf (2005) in a vegetated shear zone; while the value is about 1.0 simulated by Kang and Choi (2009) in a non-shear zone. For Run 2, the typical value of the diffusion coefficient Dt calculated by the equation (6) is 1.294E-5 m2 s−1 in the vegetation zone. If we take another alternate formula Dt = (υt + v )/St , the turbulent

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

93

Fig. 12. Computed contours of the Reynolds flux for Run 2: (a) - u c /UC0 and (b) - w c /UC0 (Dashed line is the vegetation edge).

Fig. 10. (a) Computed Reynolds stress profile for Run 2 and measured data by Okamoto and Nezu(2010), (b) Derivation of the Reynolds stress profiles with the vegetation density.

Fig. 13. Correlation coefficient with the variation of vegetation density.

5.3. Turbulence intensity and scalar flux distribution

Fig. 11. Distribution of the computed scalar turbulence intensity crms /C0 for Run 2 (Dashed line is the vegetation edge).

Schmidt number St will be 0.434, which is near to the value 0.47 by Ghisalberti and Nepf (2005). The value indicates that the concentration transport is more than twice as efficient as the momentum transport in the vegetated zone. This reveals that the present scalar flux modeled by the Eq. (6) including the mechanic diffusion term is reasonable in the vegetated flows. The simulated results therefore can reasonably produce the process of the concentration diffusion in the vegetated flows. In addition, the comparisons of the computed concentration profiles and the measured data are shown in Fig. 8. The good agreement between them shows that the present model is capable of simulating the solute concentration transport in the vegetated flows.

The mean flow and concentration fields were presented in the previous sections, whereas the turbulence statistics are reported here to investigate the turbulence and scalar turbulence structures in the vegetated flows. Fig. 9(a)–(c) depict all the normalized three components of the turbulence intensity, i.e., urms /U, vrms /U, wrms /U for Run 2, respectively. One can see from these figures that the local higher values of urms, vrms, wrms locate near the tip of vegetation edge and a steep gradient of the turbulence intensity distributes above and below the vegetation edge, because of the discontinuity of vegetation resistant drag between the vegetation and non-vegetation zone. These higher values of turbulence intensity enhance the scalar mixing in the mixing-layer region. Therefore, the difference in the three components of the turbulence intensity represents that the vegetated flow field is anisotropic, which have a great impact on the scalar mixing in the longitudinal, transversal and vertical directions. The computed Reynolds stress profile for Run 2 and the experimental data is compared in Fig. 10(a). The values are normalized by the friction velocity u∗ . The agreement between the computed result and the experimental data is generally acceptable. One can note from Fig. 10(b) that the peak of the Reynolds stress occurs near at the vegetation tip, and it decreases towards the bed and the water surface. Moreover, one can notice that the maximum value of the Reynolds stress increases with the increasing of the vegetation density.

94

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

Fig. 14. Snapshots of the coherent structures (Q = 0.05) with the concentration transport for Run 2: (a) t = 0.0 s, (b) t = 0.2 s, (c) t = 0.4 s and (d) t = 0.8 s.

Fig. 11 gives the distribution of the scalar turbulence intensity crms normalized by the value of source concentration for Run 2. It can be observed that the distribution of the scalar turbulence intensity crms is similar to the pattern observed for the time averaged concentration contours shown as Fig. 4(a), with a local peak value near the vegetation edge and a steep gradient near the

vegetation edge. It also should be noticed that the peak value of crms have a similar pattern with the maximum value of the mean concentration distribution. The peak value of crms decays and the concentration intensity distribution tends more uniform with the increase of the distance from the solution injected dye.

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

95

A scalar flux plays an important role in the scalar diffusion and mixing. Fig. 12 shows the contours of the normalized scalar flux -u c and -w c for Run 2. It can be seen that the values of the scalar flux −u c are positive and -w c are negative above the vegetation edge. The large peak values of -u c and -w c are distributed in the over-canopy region, which implies that the ejection motion transports the concentration into the above-canopy layer. Conversely, the values of -u c are negative within the canopy layer, indicating that the high distribution of concentration is transported downward by the sweep motion. The details about the mechanics of ejection-sweep motions will quantitatively be discussed in the following sections. The efficiency of scalar mixing can be calculated through the correlation coefficient, RCw = c w /σc σw , where σ C and σ w are the standard deviations of c and w, respectively. Fig. 13 presents the correlation coefficient with the variation of vegetation density. For the vegetated flows, the typical value of Rcw is about −0.3 in the mixing-layer zone, while for the non-vegetated flows, the typical value of Rcw is about −0.6. Comparison of the value of Rcw for the vegetated and non-vegetated flows reveals that the larger in magnitude of Rcw indicates that the mixing of the scalar concentration in the vegetated flows is much more efficient than those in the non-vegetated flows. The data of Fig. 13 shows that it is obvious that there is not monotonically linear relationship between the correlation coefficient and the vegetation density. 5.4. Coherent structures As we know, the LES not only provides the mean flow and concentration fields but also produces the instantaneous flow and concentration fields. In order to visualize the turbulence structure, the Q criterion is used to insight the instantaneous three dimensional coherent structures. The quantity Q is the second invariant 2 2 of the velocity gradient tensor, defined as Q = 0.5(  − S ), in which S and   are the symmetric and anti-symmetric part of the velocity gradient tensor, respectively. The quantity Q describes the relations between the strain rate and the rotation rate, and it was usually identified the coherent structures in the flows with vegetation. Fig. 14 presents the snapshots of coupling of the coherent structures with the concentration fields for Run 2. It is obvious that the evolution and development of the large-scale spanwise rollers can be observed clearly. Moreover, the small scale vortices always link with the organized large scale vortices. The life of the large-scale vortices is long and quasi-period; on the contrast, the life of the small scale vortices is short and random. It is the long and quasiperiod coherent structures, containing the major part of the kinetic energy, which contributed the most to the mass and scalar concentration transport between the vegetation and non-vegetation zone. These KH vortices were strongly accompanied by the ejection and sweep motions, which play a significant role in controlling the process of the concentration mixing. Meanwhile, the transportation of the instantaneous high scalar concentration clouds is most controlled by the quasi-periodic large scale structures, which is evidenced by the dynamics of the coherent structures with the instantaneous concentration fields shown as Fig. 14. The spanwise vortices containing the high concentration solute clouds travel downstream like a train of waves. The good synchronization can be calibrated by the power spectral densities of the record vertical instantaneous velocity component and concentration (Fig. 15(a) and (b). The peak of the power spectral density (PSD) of the vertical instantaneous velocity component is in the range of 1–2 Hz, which is the frequency of the large scale structures; while the range of the peak the PSD of the concentration is about 1–2 Hz. Obviously, there is well synchronous between the instantaneous high concentration clouds and the large

Fig. 15. PSD of the vertical instantaneous velocity and concentration located at x/hv = 3, y/hv = 4 and z/hv = 1: (a) and (b) for Run 2.

scale structures. This implies that the KH vortex is a major factor contributing to the instantaneous high concentration clouds transport. Additionally, we can observe that the KH vortices are hairpin vortexes for Run 2, as shown in Fig. 16(a). The hairpin vortexes with the -type which are usually observed in the wall boundary layer (Adrian et al., 2001) is also existence in the flow with vegetation. However, for the open channels with transitional vegetation (Cd α hv ∼0.1 Nepf, 2012), the hairpin vortexes at the layer between the vegetated and non-vegetated zone are not clearly detected as shown in Fig. 16(b). In order to investigate the changing of the vortexes topology with the variation of the vegetation density, the length scale of the vortexes based on the spatial correlation (Castro et al., 2006) is used. Fig. 17 presents the length scales of the vortexes with the variation of the vegetation density. Typically, the data show that, for the dense vegetation (Cd α hv >0.1 Nepf, 2012), the length scale of the vortexes is about 3hv, which is consistent with the experimental results (Castro et al., 2006). On the other hand, the length scales of the vortexes in the open channel with transitional vegetation or without are larger than those with dense vegetation. This is the reason of the length scale of the vortexes detected by the secondary current in the open channel

96

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

Fig. 16. Top view vortex of the iso-surface (Q = 0.4) and velocity vector (u,v,0) at the plane of z/hv = 1: (a) for Run 2 and (b) for Run 5.

Table 2 Definitions of the quadrants for scalar concentration. Quadrant

u i

c

Name

Q1 Q2 Q3 Q4

>0 <0 <0 >0

>0 >0 <0 <0

Ejection Inward interaction Sweep Outward interaction

rant of the Reynolds stress is defined as follows:



u w



i,HS

= limT →∞

T 0

u (t ).w (t )Q i,HS (u , w )dt

(i = 1, 2, 3 and 4 ),

(13)

where T is the time interval and Q i, HS (u , w ) is the detection funcFig. 17. Length scales of vortexes with the variation of vegetation density.

tion, the hole size parameter HS = 0 is used in this study. The fractional contribution Q i, HS to −u w from each event is



Q i,HS (u , w ) = with transitional vegetation or without, not by the KH vortices induced by the dense vegetation.

5.5. Quadrant analysis A quadrant conditional analysis, which provides the information on the motions responsible for momentum transport, for the fluctuating velocity is applied (Willmarth and Lu, 1972). In the quadrant analysis, there are four quadrants (i = 1– 4) in the u -w plane. For simplicity when discussing the quadrants, we will note them as Q1, Q2, Q3 and Q4. Quadrant Quadrant Quadrant Quadrant

1: 2: 3: 4:

u > 0 u < 0 u < 0 u > 0

and and and and

w > 0 w > 0 w < 0 w < 0

represents represents represents represents

outward interaction ejection inward interaction sweep

The value of –u w < 0 corresponds to a positive Reynolds stress, whereas –u w > 0 represents a negative Reynolds stress. The quad-

1 0



when |u w | > HS u w . otherwise

(14)

Here, Q i, HS > 0 is for the first and the third quadrant and Q i, HS < 0 is for the second and the forth quadrant. The definitions of the four quadrants for the scalar transport, which is widely used in the case of scalar dispersion, are altered compared to the previous definitions for momentum flux and are briefly described in Table 2. Fig. 18 presents the histogram of Q1, Q2, Q3 and Q4 for u and w with the variation of vegetation density located at x/hv = 3, y/hv = 4 and z/hv = 1. For the dense vegetation, the fraction of the ejection and sweeping events are larger than 50% while those of the inward and outward events are small than 50%, which implicates that the process is dominated by the ejection and sweeping events. The fraction in the quadrant 2 and quadrant 4 exceed those in the quadrants 1 and 3 (Fig. 19(a)), indicating the relative contributions of the ejection and sweep is dominant in the momentum transportation. We may conclude that the coherent structures are accompanied by the ejection and sweeping motions, which are

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

Fig. 18. Histogram of Q1, Q2,Q3 and Q4 for u’ against w’ located at x/hv = 3, y/hv = 4 and z/hv = 1.

major factors contributing to the Reynolds stress transportation, in the mixing-layer zone in the open channel with dense vegetation. Conversely, in the open channel with transitional vegetation, the inward and outward events are predominating. Fig. 19(b–d) presents the ratios of the inward and outward to the ejection and sweep events for C against ui  . In the streamwise direction (Fig. 19(b), the fraction of the scalar flux with vegetation occurrence within the quadrants 2 and 4 is the inward and outward events, is larger those in the quadrants 1 (ejection events) and 3 (sweeping events). The contributions of Q1 and Q3 to the scalar flux are minor. The proportion of the ejection and sweep events is little than 50% and the proportion of the inward

97

and outward events is more than 50%, respectively. This implies that the ejection and sweep events are not major in the streamwise direction. Concerning the spanwise direction, the calculated ratios of the inward and outward events to the ejection and sweep events is roughly equal to 50% (Fig. 19(c)), indicating a zero scalar flux in the lateral direction. It is a fact that the inject point is located at the central lateral plane in the open channel, indicating that the intensities of vrms and Crms are relatively uniform, respectively. Considering the vertical direction, Q1 indicates the ejection of concentration scalar clouds vertically upward away from the vegetation. In contrast, the sweeping events, as indicated by Q3, correspond to the concentration scalar clouds downward toward the bed. Ejection-sweeping phenomenon, as evidenced by the ratios of the inward and outward to the ejection and sweep events (Fig. 19(d)), is a major factor in the contribution to the scalar flux in the open channel with vegetation. They enhance the scalar mixing and create a more uniform concentration distribution in the vertical direction. However, the contributions of Q2 and Q4 to the scalar flux are minor. It is an evidence of ejection-sweeping motions, which are major contributions to the scalar flux in the mixing-layer zone. They are consistent with the fact that the ejection and sweeping are principal motions in the momentum transportation. 5.6. Spectral analysis As turbulent length scales play an important role in the eddy cascade process in the open channel flows with vegetation, the PSDs of the velocity and concentration are investigated. The computed power spectral densities of the vertical instantaneous velocity and concentration located at x/hv = 3, y/hv = 4 and

Fig. 19. Histogram of (Q2+ Q4)/(Q1 +Q3) located at x/hv = 3, y/hv = 4 and z/hv = 1: (a) u against w ,(b) c against u , (c) c against v and (d) c against w .

98

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99

6. Conclusions A 3D large eddy simulation model has been applied to simulate the flow and concentration transport in an open channel with submerged vegetation. The simulated velocity, Reynolds stress and concentration profiles are in good agreement with the experimental results by Okamoto and Nezu (2010). The typical flow and concentration distributions with vegetation can be produced by the large eddy simulation model. The velocities monotonically decrease within the vegetation and increase above the vegetation with the increasing of the plant density. The concentration monotonically decays along the vegetation edge with the increasing vegetation density. However, for the vegetation density α ≥ 30 m−1 , both the velocity changing within or above the vegetation and the scalar concentration decaying along the top of vegetation are less sensitive. The synchronization of the instantaneous high concentration clouds and the large scale structures is achieved. Moreover, the scalar transport mechanism is displayed by the dynamic of coupling of the coherent structures with the concentration fields, and is quantitatively explained by the quadrant conditional analysis. The ejection-sweeping events are dominant in the process of the momentum and scalar flux transport. Comparison of the power spectral densities computed from the time-recording velocity with vegetation or without shows that the frequency of KH vortexes generated by the vegetation is in the range of 1–2 Hz. Acknowledgments The study is supported by the National Natural Science Foundation of China (51509075, 51279047), the Post Doctoral Fund of Jiangsu Province (1501042A), the National Key Basic Research Program of China (973 Program) (2012CB417006), and the Program for Changjiang Scholars and Innovative Research Team in University (IRT1233). The authors are thankful the comments from anonymous reviewers. References

Fig. 20. PSD at point of x/hv = 3, y/hv = 4 and z/hv = 1: (a) variation of PSD of the vertical instantaneous velocity with vegetation density and (b) PSD of the instantaneous concentration for Run 1.

z/hv = 1are shown in Fig. 20. It is obvious that the PSDs of the velocity follow Kolmogorov’s −5/3 scaling for the inertial sub-range quite well. We argue that the low-frequencies (large scale structures) are corresponds to the velocity fluctuations caused by the KH vortices occurring at the interface between the non-vegetated and vegetated zones; on the other hand, the high-frequencies (small scale structures) are associated with the much smaller turbulence scale(e.g. eddies at Kolmogorov length scale). It is identified from Fig. 20(a) that the frequency of the shear layer vortices induced by vegetation is about 1.6 Hz, which matches the typical values of the KH frequency within the literature (Ghisalberti and Nepf, 2002; Siniscalchi et al., 2012). Moreover, we can find that the increasing the vegetation density has no effects on the peak of PSD. However, no obvious peak value of the PSD computed from the instantaneous velocity and concentration for Run 1 is observed as shown in Fig. 20(a) and (b). This indicates that the vegetation leads to the KH vortices instability in the vegetated shear layer zone.

Adrian, R.J., Balachandar, S., Lin, Z.C., 2001. Spanwise growth of vortex structure in wall turbulence. KSME Int. J. 15 (12), 1741–1749. Busari, A.O., Li, C.W., 2015. A hydraulic roughness model for submerged flexible vegetation with uncertainty estimation. J. Hydro Env. Res. 9 (2), 268–280. Carollo, F.G., Ferro, V., Termini, D., 2002. Flow velocity measurements in vegetated channels. J. Hydraul. Eng. 128 (7), 664–673. Castro, I.P., Cheng, H., Reynolds, R., 2006. Turbulence over urban-type roughness: deductions from wind-tunnel measurements. Boundary Layer Meteorol. 118 (1), 109–131. Choi, S.U., Kang, H.S., 2004. Reynolds stress modeling of vegetated open-channel flows. J. Hydraul. Res. 42 (1), 3–11. Cui, J., Neary, V.S., 2008. LES study of turbulent flows with submerged vegetation. J. Hydraul. Res. 46 (3), 307–316. Ghisalberti, M., Nepf, H.M., 2002. Mixing layers and coherent structures in vegetated aquatic flows. J. Geophys. Res. 107 (C2). Ghisalberti, M., Nepf, H., 2005. Mass transport in vegetated shear flows. Env. Fluid Mech. 5 (6), 527–551. Huai, W.X., Zeng, Y.H., Xu, Z.G., et al., 2009. Three-layer model for vertical velocity distribution in open channel flow with submerged rigid vegetation. Adv. Water Resour. 32 (4), 487–492. Huthoff, F., Augustijn, D., Hulscher, S.J.M.H, 2007. Analytical solution of the depth-averaged flow velocity in case of submerged rigid cylindrical vegetation. Water Resour. Res. 43 (6). Jarrin, N., Benhamadouche, S., Laurence, D., et al., 2006. A synthetic-eddy-method for generating inflow conditions for large-eddy simulations. Int. J. Heat Fluid Flow 27 (4), 585–593. Järvelä, J., 2002. Flow resistance of flexible and stiff vegetation: a flume study with natural plants. J. Hydrol. 269 (1), 44–54. Kang, H., Choi, S.U., 2009. Turbulence Modeling of Solute Transport in Open-Channel Flows Over Submerged Vegetation//Advances in Water Resources and Hydraulic Engineering. Springer Berlin Heidelberg, pp. 674–679. Leu, J.M., Chan, H.C., Jia, Y.F., He, Z.G., Wang, S.S.Y., 2008. Cutting management of riparian vegetation by using hydrodynamic model simulations. Adv. Water Resour. 31 (10), 1299–1308. Li, C.W., Xie, J.F., 2011. Numerical modeling of free surface flow over submerged and highly flexible vegetation. Adv. Water Resour. 34 (4), 468–477.

J. Lu, H. Dai / Advances in Water Resources 97 (2016) 87–99 Li, C.W., Yan, K., 2007. Numerical investigation of wave–current–vegetation interaction. J. Hydraul. Eng. 133 (7), 794–803. Lopez, F., Garcia, M., 1998. Open channel flow through simulated vegetation: suspended sediment transport modeling. Water Resour. Res. 34 (9), 2341–2352. Lu, J., Dai, H.C., 2016a. Large eddy simulation of flow and mass exchange in an embayment with or without vegetation. Appl. Math. Model. 40, 7751–7767. http://dx.doi.org/10.1016/j.apm.2016.03.026. Lu J., Dai H.C. Three dimensional numerical modeling of flow and scalar transport in vegetated channel. submit to J. Hydro Env. Res. 2016b. Lu, J., Wang, L.L., Dai, H.C., 2010. Numerical investigating of vertical turbulent jets in different types of waves. China Ocean Eng. 24 (4), 611–626. Naot, D., Nezu, I., Nakagawa, H., 1996. Hydrodynamic behavior of partly vegetated open channels. J. Hydraul. Eng. ASCE 122, 625–633. Neary, V.S., 2003. Numerical solution of fully developed flow with vegetative resistance. J. Eng. Mech. 129 (5), 558–563. Nepf, H.M., Sullivan, J.A., Zavistoski, R.A., 1997. A model for diffusion within emergent vegetation. Limnol. Oceanogr. 42 (8), 1735–1745. Nepf, H.M., Vivoni, E.R., 20 0 0. Flow structure in depth-limited, vegetated flow. J. Geophys. Res. 105 (C12), 28547–28557. Nepf, H.M.., 2012. Flow and transport in regions with aquatic vegetation. Annu. Rev. Fluid Mech. 44, 123–142. Nepf, H., Ghisalberti, M., 2008. Flow and transport in channels with submerged vegetation. Acta Geophys. 6 (3), 753–777. Nepf, H.M., 1999. Drag, turbulence, and diffusion in flow through emergent vegetation. Water Resour. Res. 35 (2), 479–489. Nezu, I., Sanjou, M., 2008. Turburence structure and coherent motion in vegetated canopy open-channel flows. J. Hydro Env. Res. 2 (2), 62–90. Okamoto, T., Nezu, I., 2010. Large eddy simulation of 3-D flow structure and mass transport in open-channel flows with submerged vegetations. J. Hydro Env. Res. 4 (3), 185–197.

99

Poggi, D., Katul, G., Albertson, J., 2006. Scalar dispersion within a model canopy: measurements and three-dimensional Lagrangian models. Adv. Water Resour. 29 (2), 326–335. Poggi, D., Porporato, A., Ridolfi, L., et al., 2004. The effect of vegetation density on canopy sub-layer turbulence. Boundary Layer Meteorol. 111 (3), 565–587. Raupach, M.R., Finnigan, J.J., Brunei, Y., 1996. Coherent eddies and turbulence in vegetation canopies: the mixing-layer analogy. Boundary Layer Meteorol. 78 (3–4), 351–382. Schultz, M.G., 2002. On the use of ATSR fire count data to estimate the seasonal and interannual variability of vegetation fire emissions. Atmos. Chem. Phys. 2 (5), 387–395. Shimizu, Y., Tsujimoto, T., 1994. Numerical analysis of turbulent open-channel flow over vegetation layer using a k–e turbulence model. J. Hydrosci. Hydraul. Eng. JSCE 11 (2), 57–67. Siniscalchi, F., Nikora, V.I., Aberle, J., 2012. Plant patch hydrodynamics in streams: mean flow, turbulence, and drag forces. Water Resour. Res. 48 (1). Spalart, P.R., Jou, W.H., Strelets, M., et al., 1997. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Adv. DNS/LES 1, 4–8. Su, X., Li, C.W., 2002. Large eddy simulation of free surface turbulent flow in partly vegetated open channels. Int. J. Numer. Meth. Fluids. 39 (10), 919–937. Tsujimoto, T., 1999. Fluvial processes in streams with vegetation. J. Hydraul. Res. 37 (6), 89–803. Werner, H., Wengle, H., 1993. Large-Eddy Simulation of Turbulent Flow Over and Around a Cube in a Plate Channel//Turbulent Shear Flows 8. Springer Berlin Heidelberg, pp. 155–168. Willmarth, W.W., Lu, S.S., 1972. Structure of the Reynolds stress near the wall. J. Fluid Mech. 55 (01), 65–92. Zhang, M.L., Li, C.W., Shen, Y.M., 2010. A 3D non-linear k–ε turbulent model for prediction of flow and mass transport in channel with vegetation. Appl. Math. Model. 34 (4), 021–1031.