Modeling of Flow Around Vegetated Meandering River Reaches

Models required in a large eddy simulation of natural river flows which are influenced by the complex bed geometry and the vegetation are proposed. First the model and the calculation method have been applied to a laboratory flow in a curved channel with artificial vegetation. The effects of vegetation on the secondary flow are confirmed to be reproduced correctly. The method is then used to simulate the flow around a real meandering river with vegetation along the banks and on the bars. The results appear reasonable and important characteristics around the bends such as varying flow velocity and the bed shear are also reproduced.


Introduction
Lower reaches of major rivers in flat land of Malaysia flow around many bends. The bed and the banks at many places are thickly covered with bushes and vegetation [1]. These two factors of river bends and the vegetation interact with each other by the mechanical processes of erosion and scouring and the vegetation succession of riparian and bar vegetation and dry land trees. The characteristics of the flow and the bed forces play the key mechanism and understanding of them is an important part of controlling and keeping the river morphology and to prevent unforeseen disasters like sudden bank failure during high flows.
The overall qualitative effects may be made by such bulk-flow parameters as the flow rate and the water stage distributions associated with possible hydrographs resulting from discharges from upstream reaches and local runoffs. However, more detailed quantities such as the cross-flow distribution and variations of flow quantities and associated sediment transport characteristics that are important to evaluate the state of the river morphology and the riparian ecosystem ( [2]- [3]). These detailed quantities can now be obtained or estimated by simulations of the detailed flow (e.g. [4]). Although one-dimensional and horizontal and vertical two-dimensional flow analyses are becoming more elaborate reflecting aspects of three-dimensional nature of the turbulent flow, the basic equations used in them do not allow inclusion of complicating effects such as those induced by strongly curving flow obstructed by submerged or emergent vegetation. The large-eddy simulation (LES), on the other hand, can simulate the details of these complex flows. Some attempts have already be done in simpler river flow situations but those considering the detailed mechanisms of the curving flows in vegetated channel still do not exist and development efforts are made.
In the present work, an LES method that the present authors have developed (Nakayama & Asami [5]) is verified against basic laboratory model flow in a vegetated curved channel. The method then is applied to simulate real flow in meandering river with thick vegetation. The bed is assumed fixed though its roughness is modelled along with the resistance die to submerged vegetation.

Basic equations and method of solution
In an application to flows in natural streams and those in coastal areas, we observe following points. First, the flow regions are very irregular in shape and change in time. Second, the scales of the flow, both spatial and timewise are very large and the spacings of grid points and the time step that can be accommodated by a readily available computer system are nowhere enough to resolve dissipation scales. Yet with a reasonable number of grid points like a few million points, which roughly divide the total flow region into cubes of orders-of-magnitude smaller in length, most of the large-scale flow and a good fraction of the turbulent fluctuations can be resolved from which not only the mean flow but turbulent stresses can be evaluated. That means that the filtered flow with this scale does represent the overall characteristics of the flow being simulated very well.
The filtered equations of motion and the continuity equation for the spatially-filtered velocity components (u,v,w) and the filtered pressure p in rectangular coordinates (x,y,z) with z taken positive vertically upward are where t is time,  and  are the density and the kinematic viscosity of water, g is the gravitational acceleration, and xx etc. are the effective stress components representing the effects of the unresolved scales of motion.
Since the movement of the water surface is an important part of the flow feature, we solve the spatially-averaged kinematic equation for the elevation h of the filtered water surface along with the velocity field. If we limit to cases where h is a single-valued function of the horizontal position (x,y), it is where us ,vs, ws are the (filtered) velocity components at the water surface, and hx and hy are the sub-grid scale free surface fluctuation terms that may not be known well but are discussed for example by Shen & Yue [6].

Sub-grid model
While a wide range of models for the sub-grid effects have been proposed in the last decade or so (e.g. Sagaut [7]), we use a model that is mathematically simple, numerically efficient and less complex in implementation is best for simulating large-scale flows in complex domains than those elaborate ones with extra computational loads. We use the eddy viscosity model with the effective kinematic viscosity sgs sgs is evaluated by the standard Smagorinsky model where  is the geometric average of the grid spacings and Cs is the Smagorinsky constant for which we typically take a value of 0.13. The numerical grid we will use is not assumed to resolve the flow near solid surfaces accurately. It is the main feature of flows of civil engineering applications that so-called near-wall flow cannot and will not be able to be resolved by any realistic computer systems that may become available in near future. How the flow in this wall region is treated is as important as how the main equations of motion (1) are solved and how the sub-grid scale effects are modeled. So we describe how the near-wall flow is treated here together with the sub-grid model.
As to the sub-grid scale terms hx and hx in the equation for the free-surface elevation, we use a gradient model ( where sgs is a model coefficient which is taken to be proportional to sgs (Hodges and Street [9]).

Model for vegetation
The vegetation and other obstacles in a stream have been modeled in various scales and various ways (e.g. Nepf [10]). The simplest is to include is as the overall resistance such as in the resistance coefficient. In LES also, the introduction of resistance is the easiest way but it can be done as three-dimensional distribution thus discriminating between the bottom roughness and the resistance due to the in-flow vegetation. When the exact shape is very complex but the envelope of the distribution may be simple enough to be resolved like trees and vegetation, the effects are modeled by introducing in-flow drag. In such a case we add to the right hand side of Eqs.(1) the body force of the following form where Cv is the average drag coefficient of the vegetation,  is the average density of the surface areas on the vegetation branches and leaves on which the drag force is exerted and V is the magnitude of the velocity at a point in the roughness influence area.

Boundary conditions
It is assumed that the flow is sub-critical at the upstream inflow and at the downstream outflow sections or take the calculation region to meet this criterion. At the inflow section, a discharge hydrograph is prescribed. First, the standard logarithmic profile for fullydeveloped two-dimensional open-channel flow is assumed at each vertical section across the width of the flow with constant free surface. The flow is then computed for certain period of time. Then the computed velocity and pressure distribution along with the free surface levels at the end of the inflow development section are fed back to the inflow section. This 'recycling' procedure allows turbulent fluctuations to develop.
In order to adjust the discharge to the given hydrograph, Q(t) the velocity component u normal to the inflow by solving equation where Qi(t) is the given discharge at time t, Si is the area of the inflow section and Ti is the relaxation time constant, which is typically set to 50 times computational time step. It is noted that when the time constant Ti is zero at all times. The free-surface elevation at the inflow section is not prescribed but is computed in the solution process.
At the outflow section, if the stage-discharge relation is known from the gauging station data, the following depth-discharge equation is used.
Here H is the mean water surface elevation at the deepest point in the outflow section and is assumed constant across the section, Qo is the total flow rate through the outflow section, CQ, r and Ho are constants taken to fit the existing data. When there is no known stagedischarge relation as in the case of the present simulation of the real river, the downstream outflow section is chosen at a position where the flow is nearly straight and nearly uniform. If the flow there is subcritical, the normal depth corresponding to the average velocity at each section in the outflow plane is given. If the flow is supercritical, the free outflow condition is applied. In the following validation calculation, the condition of subcritical uniform-flow depth is used.
The wetting edge is allowed to advances or retreats as the flow goes over dry land. No specific condition is needed here but the calculation of the free surface position should consider the possibility that the zero flow depth next to region of nonzero depth can have finite depth and vice versa.

Validation by calculation of flow in vegetated channel bend
For the validation of the present model, flow around a 60 degee bend with vegetated inner bank measured by Tominaga et al. [11] as shown in Fig.1 is calculated. This flow has also been calculated using a RANS model by Lee and Choi [12]. The total length of the channel is 17.2m consisting of 10.8m and 3.6m long straight channels upstream and downstream of the 60 degree bend. The radius of curvature of the centreline of the curved part of the channel is 2.7m and the channel width is 0.9m. The inner 1/4th of the channel width is covered by emergent vegetation. The flow condition and the vegetation density are given in Table 1. The flow around a bend is influenced by the secondary flows due to the extra acceleration of the curved flow. The correct reproduction of the secondary flows is an important part of the simulation. In this experiment, the bed slope is not given but the entire flume is adjusted to obtain steady flow with the depth of 0.15m. In the calculation an assumed bed slope is given to obtain similar velocity and the depth. Two cases with and without the vegetation zone corresponding to the experiment are simulated. Fig.2 and Fig.3 show  vegetated is towards the inner bank and is not in the same direction as the main secondary flow according to the experiment. This feature is reproduced correctly by the present simulation while the RANS results indicate the flow in the vegetated region is in the same direction as the main secondary flow. Fig.4 shows the distribution of the surface velocity at one instant. It also shows the water surface elevation distribution in color. The surface rise upstream of the vegetated region is reproduced well. The degree of velocity fluctuation is also seen reasonable.

Simulation of flow around Sungai Perak near Teluk Intan bend
Calculation was conducted for meandering river flow in a low land in Perak, Malaysia.  Some results of the present LES are shown in Fig.6 and Fig.7. Fig.6 shows the velocity vectors on the free surface and the elevation of the water surface in colored contours for the case of constant flow rate of 1000m 3 /s. The average slope of the river surface in the area is about 1/2000. The left figure is obtained without using the vegetation model and the right one is with the vegetation as shown in Fig.5. The differences are very small but the surface velocity near the vegetated banks is slightly smaller than the case without vegetation.
The main flow and the secondary flow velocity at a cross section just downstream of the first bend are shown in Fig.7. This is the area thickly vegetated and downstream of the bar. The secondary flow going towards the outer bank (left bank) is seen but near the right bank there is another flow going towards the right bank. This is downstream of the vegetated inner bank and the island and the flow is complicated. The LES calculation gives these complex secondary flows and the detailed bed shear stress though it is not shown here.

Conclusion
The resistance distributed within a flow that models the effects of vegetation has been implemented in a large eddy simulation method of flow of natural river reach. The validation calculation done for a laboratory flow in a curved channel with immergent vegetation in the inner bank shows that the overall main flow and the secondary flows are simulated well. In addition to the main secondary flow due to the bend, the outer cell which a RANS method could not reproduce is reproduced in the present LES simulation. The water level rise upstream of the vegetated zone and the restricted region of main secondary flow within the vegetation zone are also reproduced. The calculation done for a limited reach of Sungai Perak appears reasonable. The latter simulation is still on-going and will be compared with the onsite data.