Flow and transport simulation of Madeira River using three depth-averaged two-equation turbulence closure models

Li-ren YU*

1. Environmental Software and Digital Visualization (ESDV), S?o Carlos 13561-120, Brazil

2. Association of United Schools-Higher Education Center of S?o Carlos (ASSER-CESUSC), S?o Carlos 13574-380, Brazil

Flow and transport simulation of Madeira River using three depth-averaged two-equation turbulence closure models

Li-ren YU*1,2

1. Environmental Software and Digital Visualization (ESDV), S?o Carlos 13561-120, Brazil

2. Association of United Schools-Higher Education Center of S?o Carlos (ASSER-CESUSC), S?o Carlos 13574-380, Brazil

This paper describes a numerical simulation in the Amazon water system, aiming to develop a quasi-three-dimensional numerical tool for refined modeling of turbulent flow and passive transport of mass in natural waters. Three depth-averaged two-equation turbulence closure models,,, and, were used to close the non-simplified quasi-three-dimensional hydrodynamic fundamental governing equations. The discretized equations were solved with the advanced multi-grid iterative method using non-orthogonal body-fitted coarse and fine grids with collocated variable arrangement. Except for steady flow computation, the processes of contaminant inpouring and plume development at the beginning of discharge, caused by a side-discharge of a tributary, have also been numerically investigated. The three depth-averaged two-equation closure models are all suitable for modeling strong mixing turbulence. The newly established turbulence models such as themodel, with a higher order of magnitude of the turbulence parameter, provide a possibility for improving computational precision.

river modeling; numerical modeling; contaminant transport; depth-averaged turbulence models; multi-grid iterative method

1 Introduction

Almost all flows in natural rivers are turbulent. Dealing with the problems of turbulence closely related to river pollutions is challenging both for scientists and engineers, because of their damaging effects on our limited water resources and fragile environment (Yu and Salvador 2005a). It is important to develop adequate mathematical models, turbulence models, numerical methods, and corresponding analytical tools for timely simulation and predicting transport behaviors in natural and artificial waterways.

Although the significance of modeling turbulent flow and transport phenomena with high precision is clear, the numerical simulation and prediction for natural waters with complex geometry and variable bottom topography are still unsatisfying. This is mainly due to theinherent complexity of the problems that need to be considered. Any computation and simulation of flow and transport processes depend critically on the following four elements: generating a suitable computational domain with the ability to treat irregular geometrical boundaries, such as riverbanks and island boundaries; establishing applicable turbulence models; adopting efficient algorithms; and developing corresponding numerical tools (Yu and Righetto 1999).

Quasi-three-dimensional hydrodynamic models are frequently used for modeling in shallow and well-mixed waters. However, many models used in practice merely consider the turbulent viscosity and diffusivity as constants or through simple phenomenological algebraic formulas (Choi and Takashi 2000; Lunis et al. 2004; Vasquez 2005; Kwan 2009), which are to a great degree estimated according to modelers’ experience. Although other practical quasithree-dimensional hydrodynamic models are really closed by the depth-averaged two-equation turbulence closure model, they almost all concentrate on the investigations and applications of the depth-averagedmodel (Rodi et al. 1980; Chapman and Kuo 1982; Mei et al. 2002; Johnson et al. 2005; Cea et al. 2007; Hua et al. 2008; Kimura et al. 2009; Lee et al. 2011), which appeared over 30 years ago, introduced by McGuirk and Rodi in 1977. It is well known that the order of magnitude of the turbulence parameterfor themodel is very low.

Recent developments in turbulence model theory have provided more realistic and advanced models of turbulent flow. From an engineering perspective, two-equation turbulence closure models can establish a higher standard for numerical approximation of main flow behaviors and transport phenomena in terms of efficiency, extensibility, and robustness (Yu and Righetto 1998, 2001). However, the most common standard two-equation closure models, used widely by various industrial departments, cannot be directly employed in quasithree-dimensional modeling. The depth-averaged versions of corresponding turbulence models should be established in advance.

Except for the closure of the classical depth-averagedmodel, current simulations still adopt the closures of the depth-averagedmodel and of depth-averagedmodel. The depth-averagedmodel stemmed from the most common standard -kωmodel, originally introduced by Saffman but popularized by Wilcox (1998). The results, computed by the three depth-averaged two-equation turbulence closure models, are compared with each other in the paper. Such example, however, hardly exists for the simulation of contaminant transport in natural waters. Modeling using different two-equation closure approaches certainly increases the credibility of simulation results (Yu and Yu 2009a, 2009b).

On the other hand, recent advancements in grid generation techniques, numerical methods, and information technology have provided suitable approaches to generate non-orthogonal boundary-fitted coordinates with collocated grid arrangement, with which the non-simplified hydrodynamic fundamental governing equations can be solved by themulti-grid iterative method (Ferziger and Peric 2002). This paper describes a quasi-threedimensional hydrodynamic simulation of flow and transport in a rather complex domain, which aimed to develop thegrid-generator,flow-solver, andgraphical user interface(orinterface), and finally formed a software. The developed software, named Q3drm1.0 by the author, provides three selectable depth-averaged two-equation turbulence closure models and can solve quasi-three-dimensional refined flow and transport phenomena in various complex natural and artificial waterways.

2 Fundamental hydrodynamic governing equations

The complete, non-simplified fundamental governing equations of quasi-threedimensional computation, in terms of coordinate-free vector forms derived using the vertical Leibniz integration method for a control volume (CV, an arbitrary quadrilateral with a center pointP), considering the variation of the bottom topography and water surface and neglecting minor terms in the depth-averaging procedure, can be written as follows:

wheretis time,Ωis the CV’s volume,Sis the face,vis the depth-averaged velocity vector,nis the unit vector in the normal direction to the face, the superscript “?” indicates that the value is strictly depth-averaged,is any depth-averaged conserved intensive property (for mass conservation,; for momentum conservation,is the components in different directions ofv; and for the conservation of a scalar,is the conserved property per unit mass),Γis the diffusivity fordenotes the source or sink of,his the local water depth atP, andρis the density.

For the momentum conservation of Eq. (1),(the depth-averaged effective viscosity); for temperature or concentration transport,(temperature or the concentration diffusivity), where the superscript “～” indicates the quantity characterizing depth-averaged turbulence. The source or sink termfor momentum conservation may include surface wind shear stresses, bottom shear stresses, pressure terms, and additional point sources (or point sinks).

The continuity and momentum equations as well as the transport equation of the scalarhave been reported in detail by Yu and Yu (2009a).

3 Depth-averaged turbulence closure models

The depth-averaged effective viscosityand diffusivity, appearing in Eq. (1), are dependent on the molecular dynamic viscosityμand depth-averaged eddy viscosity:and, whereσφ,tis the turbulence Prandtl number for temperature diffusion or the Schmidt number for concentration diffusion, andis a scalar, usually determined by two extra transported variables.

The first two-equation turbulence model for depth-averaged calculation was suggested byMcGuirk and Rodi (1977):

whereandstand for the depth-averaged turbulent kinetic energy parameter and the dissipation rate parameter of, respectively,andare the source or sink terms,Pkis the production of turbulent kinetic energy due to interactions of turbulent stresses with horizontal mean velocity gradients, and the turbulent viscositycan be expressed as

The values of empirical constantsCμ,σk,σε,C1, andC2in Eqs. (2) through (4) are the same as in the standardk-εmodel, equal to 0.09, 1.0, 1.3, 1.44, and 1.92, respectively. The additional source termsPkvandPεvin Eqs. (2) and (3), respectively, are mainly produced by the vertical velocity gradients near the bottom, and can be expressed as follows:

where the local friction velocityu*is equal to, whereandare the depth-averaged velocity components in thexandydirections, respectively;Cfrepresents an empirical friction factor; and the empirical constantsCkandCεfor open channels and rivers are

wheree*is the dimensionless diffusivity of the empirical formula for undisturbed channel and river flowswithU*being the global friction velocity.

In 1989, Yu and Zhang developed a depth-averaged second-order closure model,, which is originated from the revisedk-wmodel developed by Ilegbusi and Spalding (1982) and was adopted as the second turbulence closure model in this study. The turbulence parameter equations (equation andequation) are

whereis the depth-averaged time-mean-square vorticity fluctuation parameter of turbulence;is the source or sink term;whereLis the characteristic distance of turbulence andxiis the coordinate (i= 1, 2); andΩstands for the mean movement vorticity. In themodel, the turbulent viscosity is defined as

The equations of turbulence parameters (equation andequation) should be solved in this model as well. The values of empirical constantsCμ,σk,σw,C1w,C2w,, andC3ware the same as those of the standardk-wmodel, equal to 0.09, 1.0, 1.0, 3.5, 0.17, 17.47, and 1.12, respectively. The corresponding additional source termsPkvandPwv, also mainly produced by the vertical velocity gradients near the bottom, can be expressed as

which were provided by Yu and Zhang (1989). The empirical constantsCwfor open channels and rivers can be written as

Recently, the author has established a new depth-averaged model,, based on the most common standard -kωmodel (ωis the special dissipation rate). The standard -kωturbulence model has been used in engineering research (Riasi et al. 2009; Kirkgoz et al. 2009). In the depth-averagedturbulence model, the turbulent viscosity is expressed as

whereis the special dissipation rate parameter of turbulence kinetic energy in the depth-averaged sense. As the third used two-equation turbulence closure model in the study,andare determined also by solving two turbulence parameter transport equations, theequation andequation (Yu and Yu 2009a):

while the empirical constantCωfor open channels and rivers can be expressed as

Except for the widely investigated and applied classical depth-averagedturbulence model, the author also adopts theandmodels mentioned above to close the fundamental governing equations in this study. The mathematical model and turbulence models developed by the author have been numerically investigated with laboratory and site data for different flow situations (Yu et al. 1990; Yu and Righetto 1998, 1999, 2001; Yu and Salvador 2005b; Yu et al. 2007; Yu and Yu 2009a, 2009b; Yu 2009). In the established mathematical models, the originally suggested empirical constants of three depth-averagedturbulence models are employed and have never been changed.

4 Grid generation

One reach of the Madeira River in Brazil, near Novo Aripuan? City situated at the middle reach of the river, has been studied by the developedgrid-generator, written in the Fortran language. In this simulation, one tributary feeds into the river reach on the right side. The confluent tributary has a concentration difference in comparison with the mainstream, caused by the humus in tropical rain forests (produced by tropic rains). With the help of the developedinterface, it is possible to determine the scale of the digital map (the Google map), to conveniently collect geometrical data, including the positions of two curved riverbanks and two boundaries of one island as well as the location of the tributary in the computational domain, and finally to generate one text file. This file contains all messages, which illustrate necessary control variables and characteristic parameters, including those on four exterior boundaries (the southeast inlet section, northeast outlet section, and the south and north riverbanks) and two interior boundaries (south and north boundaries of the island), and can be read bygrid-generatorfor generating the expectant coarse and fine grids (two-level grids).

Fig. 1 displays the digital map, on which the developedinterfacehas collected the location messages of the riverbanks and the island and divided the computational river reach into 47 sub-reaches with 48 short cross-river lines. It should be noted that the cross-river lines between the riverbank and island boundary have been redrawn, in order to involve the island configuration. Fig. 2(a) presents the generated non-orthogonal body-fitted coarse grid, with a resolution of 228 nodal points in theidirection and 18 nodal points in thejdirection, respectively. In the generated grid, the nodal points in the transversal grid lines are uniform. The total length of the calculated river reach is 19.606 km.The flow direction is from the southwest to the northeast. The tributary feeds into the mainstream on the south riverbank, with the numbers of nodal points fromi= 66 toi= 74 on the coarse grid. The island starts at (i= 56,j= 9) and ends at (i= 118,j= 9) on the same grid. The developedgrid-generatorgenerated two layers’ grids, on which all of the geometric data necessary in the later calculation of flow and transport must be stored and then can be read by the developedflow-solver. The resolution of the fine grid is 454 × 34, and it is displayed in Fig. 2(b). This means that one grid cell of the coarse grid was divided into four grid cells of the fine grid.

5 Solutions of flow and side discharge

The behaviors of flow and contaminant transport were simulated with the developedflow-solver(also written in the Fortran language), in which the SIMPLE (semi-implicit method for pressure-linked equation) algorithm for FVA (finite volume approach), Guass’divergence theorem, ILU (incomplete lower-upper) decomposition, PWIM (pressure weighting interpolation method), SIP (strongly implicit procedure), under relaxation, and multi-grid iterative method have been used. The fundamental governing equations were solvedat each grid level with the following: two momentum equations (theequation andequation), one pressure-correction equation (theequation), one concentration transport equation (theequation, whereis the concentration), and two transport equations of turbulence parameters (theequation andequation; orequation andequation; orequation andequation).

Fig. 1 Digital map of study reach

Fig. 2 Generated computational grids

The calculated main stream flow rate in low water season is 3500 m3/s, while the width, area, and mean water depth of the inlet section are 732.7 m, 3571.9 m2, and 4.88 m, respectively.The empirical friction factor (Cf) equals 0.004 98. The tributary flow rate and the concentration difference of humus (ΔC) are 100 m3/s and 50 mg/L, respectively. The velocity components and turbulence parameters at the main stream inlet and tributary inpouring section are uniform. On the outlet section, the variables satisfy the constant gradient condition. Three depth-averaged two-equation turbulence closure models,, andmodels, were adopted to close the quasi-three-dimensional hydrodynamic model. Thecorresponding turbulence parameters of three turbulence models can be calculated by empirical formulas:, andat the main stream inlet sections are 0.046 5 m2/s2, 0.000 96 m2/s3, 0.257 s?2, and 0.23 s?1, respectively, and, andat the tributary inpouring section equal 0.003 29 m2/s2, 0.000 026 m2/s3, 0.015 3 s?2, and 0.087 6 s?1, respectively. The wall function approximation has been used to determine the values of velocity components and turbulence parameters at the nodal points in the vicinity of riverbanks and island boundaries.

Due to the existence of the island in the mesh, the values of the under-relaxation factors for velocity components, pressure, concentration, and two turbulence parameters in the multi-grid iterative method may be lower than those with no island existing in the domain. However, in this example, these factors are still 0.6, 0.6, 0.1, 0.7, 0.7, and 0.7, respectively. The maximum allowed numbers of inner iterations for solving velocity components, pressure, concentration, and two turbulence parameters are 1, 1, 20, 1, 1, and 1, respectively. The convergence criteria for inner iteration are 0.1, 0.1, 0.01, 0.1, 0.01, and 0.01, respectively. The value of the Stone’s solver is equal to 0.92. The normalized residuals for solving velocity, pressure, concentration, and turbulence parameter fields are all less than the pre-determined convergence criterion (1.0 × 10?4).

The simulation obtained various distributions of flow, pressure, concentration, and the turbulence parameter, which are useful to the analyses of interesting problems in engineering. Some of the results, simulated with the, andmodels on the fine grid, are presented in Figs. 3 through Fig. 7.

Fig. 3 Results calculated withmodel

Fig. 4 Three-dimensionalfields, calculated with, andmodels

Fig. 5 Three-dimensional, andfields

Fig. 6 Three-dimensionalfields, calculated with, andmodels

Fig. 7 Concentrations along south riverbank

Fig. 3 displays the results calculated with themodel. Fig. 3(d) illustrates that thepollutant plume well develops along the right riverbank at the lower reach of the inpouring section. The distributions of the same depth-averaged physical variables and turbulence parametercalculated with theandmodels are similar to those in Figs. 3(a) through (e). Figs. 4(a), (b), and (c) demonstrate the three-dimensional distributions of turbulence parametercalculated with the three depth-averaged turbulence models. They are quite similar to one another, with the following maximum values: 0.264 4 m2/s2formodeling (Fig. 4(a)), 0.264 5 m2/s2formodeling (Fig. 4(b)), and 0.269 9 m2/s2formodeling (Fig. 4(c)). Figs. 5(a), (b), and (c) present the three-dimensional distributions of the turbulence parameters, and, which are different from one another, because of the different definitions of the second turbulence parameter. The maximum value ofin Fig. 5(a) is 0.003 4 m2/s3; the same values ofin Fig. 5(b) and ofin Fig. 5(c) are 0.257 4 s?2and 0.502 s?1, respectively. Figs. 6(a), (b), and (c) illustrate the three-dimensional distributions of turbulent viscosity, calculated using Eq. (4) formodeling (Fig. 6(a)), Eq. (9) formodeling (Fig. 6(b)) and Eq. (12) formodeling (Fig. 6(c)). Basically, they are similar to one another, especially forandmodeling, where the maximum values ofare 2 803 Pa·s (Fig. 6(a)) and 2 802 Pa·s (Fig. 6(b)), respectively; but the maximum value formodeling is 2 817 Pa·s (Fig. 6(c)). Fig. 7 displays the comparisons of concentration profiles on the fine grid, with (a) being the profile along the centers of the grid cells atj= 2 (i.e., along a curved line from the inlet to the outlet near the south riverbank), (b) being the transversal concentration profile ati= 150 (i.e., downstream close to the tributary outlet), and (c) being the transversal concentration profile ati= 454 (i.e., at the main stream outlet), calculated with the depth-averaged, andturbulence models, respectively. Fig. 8 shows the comparisons of the turbulence parameters, andon the transversal section ati= 300 of the fine grid, also calculated with the three depth-averaged turbulence models.

Fig. 8 Turbulence parameters, andati= 300

6 Plume development at beginning of discharge

In order to understand the development process of a contaminant plume, a special simulation was performed using themodel closure for the following case. The concentration of the tributary first equals zero, and then the value of concentration instantaneously reaches 50 mg/L att= 0, while the flow rates, either of main stream or of tributary, remain constant. Figs. 9(a) through (h) illustrate the plume development and thevariation in the lower reach of the tributary outlet section, of which Fig. 9(a) shows the clean water confluence, and Figs. 9(b) through (h) display the process of contaminant inpouring and plume development with a constant time difference ΔT(4 500 s) from each other.

Fig. 9 Plume development

7 Discussion and conclusions

Two levels of grids were used in this simulation: a coarse grid and a fine grid. The simulation on these two grids can satisfy the simulation demand. For example, by setting the number of grid levels at three, computations not only on coarse and fine grids but also on finest grid can be realized. The selection of the number of grid levels depends on the solved problems and computational requirements.

The uniform distributions of variables at the main stream inlet and tributary inflow section were used in this simulation, considering that the distance from the inlet to the tributary is long enough and the confluent flow rate from the tributary is rather small compared with the main stream flow rate. Improvement is possible through, for example, using multiple block technology or involving a part of the tributary in the computational domain. These intentions should be tested, after obtaining the site data through an expected site observation.

The solved depth-averaged concentration variable in the current computation is the concentration difference of humus between the confluent tributary black water and the main stream clean water (50 mg/L). Some contaminant indexes of discharged black water, such as COD and BOD, can also be considered the solved variable. The developed numerical tool of this study possesses the ability to simultaneously solve two concentration components in one calculation, which are caused by industrial, domestic, and natural discharges.

Fig. 4 demonstrates that the distributions of the turbulence parameter, calculated with the three turbulence models, vary strongly in the computation domain, but are quite similar to one another. However, the characteristics of the distributions of turbulence parameters, and, shown in Figs. 5(a), (b), and (c), respectively, are different from one another, though they also vary sharply. The calculated eddy viscosity, presented in Figs. 6(a), (b), and (c), also varies strongly. In fact, the eddy viscosity changes from point to point in the computation domain, especially in the areas near riverbanks and island boundaries. To solve the problems of contaminant transport caused by side discharge, for example, the pollution plume usually develops along a region near the riverbank (Fig. 3(d)), whereactually varies much strongly (Fig. 6). This means thatshould be precisely calculated using suitable higher-order turbulence closure models with higher precision, and cannot be considered an adjustable constant.

Fig. 7(a) shows that the computational concentration profiles along the south riverbank, either from themodel, from themodel, or from themodel, only have a quite small difference from one another, in the range of 0.02 mg/L to 0.5 mg/L. This means that the three utilized depth-averaged two-equation turbulence closure models almost have the same ability to simulate contaminant plume distributions along the riverbank. This conclusion also coincides with the result of the previous research that the depth-averaged two-equation turbulence closure models are all suitable for modeling strong mixing turbulence (Yu andRighetto 2001). However, the abilities and behaviors of different depth-averaged two-equation turbulence closure models for rather weak mixing, also often encountered in engineering, should be further investigated.

Figs. 7(a), (b), and (c) also illustrate that the peak values of the depth-averaged concentration plume decreases along the longitudinal direction, i.e., from 34 mg/L ati= 150 to 7.5 mg/L ati= 454, with a minor increment of transversal width (i.e., the distance from the bank to the given concentration contour on transversal section), owing to the transversal diffusion effect, for example, 210 m for the 1.0 mg/L contour ati= 150 and 305 m for the 1.0 mg/L contour ati= 305, displayed in Figs. 7(b) and (c), respectively. It should be noted that in natural waters, not only the velocity gradient but also the variation of bottom topography and other factors may turn into the main factors affecting plume configuration. Actually, in this calculation the water depth varied quickly from 1.05 m (nearest the bank grid point) to 3.1 m (whereCΔ= 1.0 mg/L) ati= 150, and from 1.05 m to 4.0 m ati= 454. Without a doubt, the increment of water depth obviously reduces the transversal enlargement of the plume, shown in Fig. 9. This feature probably is rather different from the plume observation in the laboratory channel with the flat bottom.

Except for the different definitions of turbulence parameters, and, Fig. 8 numerically demonstrates that the order of the magnitude of the turbulence parameteris smaller than the order of the magnitude of, and much smaller than the order of the magnitude of. In the simulation, thevalue ranges only from 1.119×10?5m2/s3to 0.003 4 m2/s3. However, theparameter andparameter range from 1.72 × 10?4s?2to 0.257 s?2and from 1.29 × 10?2s?1to 0.502 s?1, respectively. Three turbulence parameters,, andall appear in the denominators of Eqs. (4), (9), and (12), which were used to calculate turbulent eddy viscosity. For numerical simulation, the occurrence of numerical error is unavoidable, especially in place near irregular boundaries. It is clear that a small numerical error, caused by solving theequation, for example, will bring larger errors for calculating eddy viscosity than the same error caused by solving the other two equations (theequation andequation) will. Without a doubt, the elevation of the order of magnitude of the second turbulence parameter, reflecting the advance of two-equation closure models, provides a possibility for users to improve their computational precision. The insufficiency of the classical depth-averagedturbulence model may be avoided by adopting other turbulence models that have appeared recently, such as themodel.

The developedGraphical User Interfaceof Q3drm1.0 software can be used in various Windows-based microcomputers. The pre- and post-processors of this numerical tool, supported by a powerful self-contained map support tool together with a detailed help system, can help the user to easily compute the flows and contaminant transport behaviors in natural waters, closed using three depth-averaged two-equation turbulence models, and to draw and analyze two- and three-dimensional graphics of computed results.

Cea, L., Pena, L., Puertas, J., Vázquez-Cendón, M. E., and Pe?a, E. 2007. Application of several depth-averaged turbulence models to simulate flow in vertical slot fishways.Journal of Hydraulic Engineering, 133(2), 160-172. [doi:10.1061/(ASCE)0733-9429(2007)133:2(160)]

Chapman, R. S., and Kuo, C. Y. 1982.A Numerical Simulation of Two-Dimensional Separated Flow in a Symmetric Open-Channel Expansion Using the Depth-Integrated Two Equation (k-ε) Turbulence Closure Model. Blacksburg: Department of Civil Engineering, Virginia Polytechnic Institute and State University.

Choi, M., and Takashi, H. 2000. A numerical simulation of lake currents and characteristics of salinity changes in the freshening process.Journal of Japan Society of Hydrology and Water Resources, 13(6), 439-452. (in Japanese)

Ferziger, J. H., and Peric, M. 2002.Computational Methods for Fluid Dynamics. 3rd ed. Berlin: Springer.

Hua, Z. L., Xing, L. H., and Gu, L. 2008. Application of a modified quick scheme to depth-averagedk-εturbulence model based on unstructured grids.Journal of Hydrodynamics, Ser. B, 20(4), 514-523.

Ilegbusi, J. O., and Spalding, D. B. 1982.Application of a New Version of the k-w Model of Turbulence to a Boundary Layer with Mass Transfer. London: Imperial College of Science and Technology.

Johnson, H. K., Karambas, T. V., Avgeris, I., Zanuttigh, B., Gonzalez-Marco, D., and Caceres, I. 2005. Modelling of waves and currents around submerged breakwaters.Coastal Engineering, 52(10-11), 949-969.

Kimura, I., Uijttewaal, W. S. J., Hosoda, T., and Ali, M. S. 2009. URANS computations of shallow grid turbulence.Journal of Hydraulic Engineering, 135(2), 118-131. [doi:10.1061/(ASCE)0733-9429(2009) 135:2(118)]

Kirkgoz, M. S., Akoz, M. S., and Oner, A. A. 2009. Numerical modeling of flow over a chute spillway.Journal of Hydraulic Research, 47(6), 790-797. [doi:10.3826/jhr.2009.3467]

Kwan, S. 2009.A Two Dimensional Hydrodynamic River Morphology and Gravel Transport Model. M. S. Dissertation. Vancouver: University of British Columbia.

Lee, J. T., Chan, H. C., Huang, C. K., Wang, Y. M., and Huang, W. C. 2011. A depth-averaged two-dimensional model for flow around permeable pile groins.International Journal of the Physical Sciences, 6(6), 1379-1387.

Lunis, M., Mamchuk, V. I., Movchan, V. T., Romanyuk, L. A., and Shkvar, E. A. 2004. Algebraic models of turbulent viscosity and heat transfer in analysis of near-wall turbulent flows.International Journal of Fluid Mechanics Research, 31(3), 60-74. [doi:10.1615/InterJFluidMechRes.v31.i3.60]

McGuirk, J. J., and Rodi, W. 1977.A Depth-Averaged Mathematical Model for Side Discharges into Open Channel Flow. Karlsruhe: Universit?t Karlsruhe.

Mei, Z., Roberts, A. J., and Li, Z. Q. 2002. Modelling the dynamics of turbulent floods.SIAM Journal on Applied Mathematics, 63(2), 423-458.

Riasi, A., Nourbakhsh, A., and Raisee, M. 2009. Unsteady turbulent pipe flow due to water hammer usingk-ω turbulence model.Journal of Hydraulic Research, 47(4), 429-437. [doi:10.1080/00221686.2009. 9522018]

Rodi, W., Pavlovic, R. N., and Srivatsa, S. K. 1980. Prediction of flow pollutant spreading in rivers.Transport Models for Inland and Coastal Waters,Proceedings of the Symposium on Predictive Ability, 63-111. Berkeley: University of California Academic Press.

Vasquez, J. A. 2005.Two Dimensional Finite Element River Morphology Model. Ph. D. Dissertation. Vancouver: University of British Columbia.

Wilcox, D. C. 1998.Turbulence Modeling for CFD. La Canada: DCW Industries, Inc.

Yu, L. R., and Zhang, S. N. 1989. A new depth-averaged two-equation (-k w~~) turbulent closure model.Journal of Hydrodynamics, Ser. B, 1(3), 47-54.

Yu, L. R., Zhang, S. N., and Cai, S. T. 1990. A new model and its numerical investigation for turbulent depth-averaged two-equation simulation in water environment.Journal of Hydraulic Engineering, 21(3), 13-21. (in Chinese)

Yu, L. R., and Righetto, A. M. 1998. Tidal and transport modeling by using turbulencek~-w~ model.Journal of Environmental Engineering, 124(3), 212-221. [doi:10.1061/(ASCE)0733-9372(1998)124:3(212)]

Yu, L. R., and Righetto, A. M. 1999. Turbulence models and applications to natural waters.Silva, R. C. V., ed.,Numerical Methods in Water Resource, Vol. 4, Chapter I, 1-122. Rio de Janeiro: Brazilian Association of Water Resources (ABRH). (in Portuguese)

Yu, L. R., and Righetto, A. M. 2001. Depth-averaged turbulencek~-w~ model and applications.Advances in Engineering Software, 32(5), 375-394. [doi:10.1016/S0965-9978(00)00100-9]

Yu, L. R., and Salvador, N. N. B. 2005a. Modeling water quality in rivers.American Journal of Applied Sciences, 2(4), 881-886. [doi:10.3844/ajassp.2005.881.886]

Yu, L. R., and Salvador, N. N. B. 2005b. Flow modelling for natural river.Proceedings of the 2nd International Yellow River Forum on Keeping Healthy Life of the River, Volume III (ISBN7-80621-999-4), 202-207. Zhengzhou: International Yellow River Forum (IYRF).

Yu, L. R., Salvador, N. N. B., and Yu, J. 2007. Models, methods and software development of contaminant transport in meandering rivers.Proceedings of the 3rd International Yellow River Forum on Sustainable Water Resources Management and Delta Ecosystem Maintenance, Volume VI (ISBN 978-7-80734–296-0), 47-56. Dongying: International Yellow River Forum (IYRF).

Yu, L. R. 2009.The First Year’s Technical Report of CFD Software Development for Refinedly Modeling Quasi-3D Flow and Contaminant Transport in Open Channels and Natural Rivers. S?o Carlos: ESDV (Environmental Software and Digital Visualization). (in Portuguese)

Yu, L. R., and Yu, J. 2009a. Numerical research on flow and thermal transport in cooling pool of electrical power station using three depth-averaged turbulence models.Water Science and Engineering, 2(3), 1-12. [doi:10.3882/j.issn.1674-2370.2009.03.001]

Yu, L. R., and Yu, J. 2009b. Quasi 3-D refined flow and contaminant transport software and its application in river water mixing in the Amazon River.Proceedings of the 4th International Yellow River Forum on Ecological Civilization and River Ethics,Volume II (ISBN 978-7-80734-797-2). Zhengzhou: International Yellow River Forum (IYRF).

This work was supported by FAPESP (Foundation for Supporting Research in S?o Paulo State), Brazil, of the PIPE Project (Grant No. 2006/56475-3).

*Corresponding author (e-mail:lirenyu@yahoo.com)

Received Jul. 25, 2011; accepted Dec. 15, 2011