ENVIRONMENTAL FLOW AND CONTAMINANT TRANSPORT MODELING IN THE AMAZONIAN WATER SYSTEM BY USING Q3DRM1.0 SOFTWARE

This paper reports a fine numerical simulation of environmental flow and contaminant transport in the Amazonian water system near the Anamã City, Brazil, solved by the Q3drm1.0 software, developed by the Authors, which can provide the different closures of three depth-integrated two-equation turbulence models. The purpose of this simulation is to refinedly debug and test the developed software, including the mathematical model, turbulence closure models, adopted algorithms, and the developed general-purpose computational codes as well as graphical user interfaces (GUI). The three turbulence models, provided by the developed software to close nonsimplified quasi three-dimensional hydrodynamic fundamental governing equations, include the traditional depth-integrated two-equation turbulence model, the depth-integrated two-equation turbulence model, developed previously by the first Author of the paper, and the depth-integrated two-equation turbulence model, developed recently by the Authors of this paper. The numerical simulation of this paper is to solve the corresponding discretized equations with collocated variable arrangement on the non-orthogonal body-fitted coarse and fine two-levels’ grids. With the help of Q3drm1.0 software, the steady environmental flows and transport behaviours have been numerically investigated carefully; and the processes of contaminant inpouring as well as plume development, caused by the side-discharge from a tributary of the south bank (the right bank of the river), were also simulated and discussed in detail. Although the three turbulent closure models, used in this calculation, are all applicable to the natural rivers with strong mixing, the comparison of the computational results by using the different turbulence closure models shows that the turbulence model with larger turbulence parameter provides the possibility for improving the accuracy of the numerical computations of practical problems.


Introduction
Numerous environmental flows can be considered as shallow, also known as quasi threedimensional flow, i.e., the horizontal length scales of the flow domain are much larger than the depth. Except for meteorological applications, typical examples are found in open channel, the middle and lower reaches of rivers, lowland rivers, lakes and reservoirs, estuaries and coastal arears as well as in oceanic flows.
The phenomena of environmental flow and pollutant transport are closely related to the human activities and social life in the densely populated and economically developed regions. People often need to build the depth-integrated (depth-averaged) mathematical models and corresponding turbulence closure models for simulating and predicting the flow and pollutant transport behaviors in shallow waters.
Almost all of flows in engineering are turbulence. Two-equation turbulence model closure can establish a high standard, which can be achieved in engineering, for the numerical approximation on main behaviors and transport processes, i.e., the so-called refined simulation. Although depthintegrated mathematical models in association with relative two-equation closure models are often used in well-mixed shallow waters, unfortunately, there still exist some models, in which the depth-integrated turbulent viscosity and diffusivity were only considered as constants or simple phenomenological algebraic formulas [1][2][3], the values of them were estimated by the modellers' experiences. On the other hand, people seem to have become accustomed to investigating and using the traditional and simple depth-integrated   k model [4][5][6][7], which has been around for 40 years, but other depth-integrated turbulence models have not been thoroughly studied and discussed. It is well known that the order of magnitude of transported variable  in   k model is very low indeed.
The turbulence modeling theory has provided quite realistic closure models, such as RANSbased models, DNS and LES, etc.
[8], however, from an engineering point of view two-equation closure turbulence models, i.e. one class of RANS-based models, can build a higher standard for numerical approximation. It is unfortunate that the so-called 'standard' two-equation turbulence model, widely adopted by various Computational Fluid Dynamics (CFD) computing engines in industry, cannot be directly adopted in depth-integrated modeling, they have to be modeled further, based on deeply understanding turbulence model theory, in order to establish the practical depth-integrated turbulence new models. The depth-integrated simulation also called depth-averaged simulation, quasi 3-D simulation, shallow water simulation as well as two and half dimensional simulation. integrated   k model was stemmed from the most common 'standard' k-ω model, originally introduced by Saffman [9] but popularized by Wilcox [10]. In this paper, the computational results, by using three turbulence model closures were carefully compared and discussed. In fact, people rarely use different depth-integrated turbulence models to explore a practical problem of the flow and transport modeling in natural shallow waters. It is clean that modeling by using different depth-integrated two-equation closure models will effectively enrich people's cognition of refined closure models, and will greatly increase the credibility of the obtained numerical results.

The depth-averaged turbulence
This paper reports a numerical simulation of environmental flow and contaminant transport in the Solimões River near The Anamã City, Brazil, aiming to develop an engineering software, namely Q3drm1.0, composed of grid-generator, flow-solver (computing engine), GUI and help system etc., which can provide three selectable depth-integrated two-equation closure turbulence models, to refinedly solve quasi 3D flow and contaminant transport problems in complex shallow waters. In the developed software, recent advancements in grid generation technologies, numerical algorithms and IT techniques have been adopted to generate non-orthogonal boundary-fitted coordinates with collocated grid arrangement, on which the fundamental governing equations can be closed by selected depth-integrated turbulence models, and solved by multi-grid iterative acceleration method [11].

Hydrodynamic Fundamental Governing Equations
Considering the variation of the bottom topography and water surface, and neglecting minor terms in the depth-integrated procedure, the complete hydrodynamic fundamental governing equations of quasi 3D computation, in terms of coordinate-free vector forms derived by using vertical Leibniz integration for a Control Volume (CV, an arbitrary quadrilateral with center point P), can be expressed as: where v is the depth-integrated velocity vector; the superscript "  " indicates that the value is strictly depth-averaged;  is any depth-integrated conserved intensive property (for mass conservation, 1   ; for momentum conservation,  is the components in different directions of v ; for conservation of a scalar,  is the conserved property per unit mass);  is the CV's volume; S is the face; q  denotes the source or sink of  ;  is the diffusivity for the quantity  ; and h and  are local water depth at P and density, respectively.

Depth-Integrated Turbulence Closure Models
In , based on the 'standard' k - model (in which ω is the special dissipation rate), originally introduced by Saffman (1970) but popularized by Wilcox (1998). The 'standard' k - turbulence model, however, has been used in engineering researches [12]. In depth-integrated   k model, the turbulent viscosity is expressed by: Where k and  stand for the depth-integrated turbulent kinetic energy and special dissipation rate of k in the depth-integrated sense. These two variables are determined by solving the k -eq. and  -eq., which can be expressed as follows [13]: is the production of turbulent kinetic energy due to the interactions of turbulent stresses with horizontal mean velocity gradients. According to the dimensional analysis, the additional source terms kv P in k-Eq. (3) and v P  in ω-Eq. (4) are mainly produced by the vertical velocity gradients near the bottom, and can then be expressed as follows: While the local friction velocity u * is equal to Where C f represents an empirical friction factor and e* is the dimensionless diffusivity of the empirical formula for undisturbed open channel/river flows  t =e*U * h with U * being the global friction velocity.
The author also uses depth-integrated  [14] as early as in 1977: Where k S and  S are the source-sink terms, t  is expressed as: Where  stands for dissipation rate of k . The values of empirical constants C  ,  k ,   , 1 C and 2 C in Eqs. (7)(8)(9) are the same as the 'standard' k-ε model, i.e. equal to 0.09, 1.0, 1.3, 1.44 and 1.92, respectively. The additional source terms P kv and P εv in Eqs. (7) And (8) can be expressed as: Where the empirical constants, C k and C ε for open channels and rivers, are: The third adopted depth-integrated second-order closure model in the paper, i.e. w k , was also developed by the author of the paper and his colleague previously [15]. The model originated from the revised kw model developed by Ilegbusi and Spalding [16]. Two extra transport equations of this model (i.e., the k -Eq. and the w -Eq.) are: and L is the characteristic distance of turbulence;  stands for mean movement vorticity; k S and w S are the source-sink terms. In w k model, the turbulent viscosity is defined as: Where w is depth-integrated time-mean-square vorticity fluctuation of turbulence. The transport equations (the k -Eq. and w -Eq.) should be solved in the model as well. The empirical constant are the same as those of 'standard' k-w model, i.e., equal to 0.09, 1.0, 1.0, 3.5, 0.17, 17.47 and 1.12, respectively. The corresponding additional source terms P kv and P wv , also mainly due to the vertical velocity gradients near the bottom, and can be written as: The empirical constant C w for open channels and rivers can be expressed as: The developed mathematical model and turbulence models in the paper have been investigated numerically by using laboratorial and site data for different flow situations [15,17]. For the Amazon River system, one example, shown in Fig.1, illustrates a comparison between the outline of black-water plume, shown on the Google satellite map, and the fine light-blue concentration contour with 35mg/L, calculated by using   k model on fine grid and plotted by the Field Browser of Q3drm1.0. In this computation, one reach of the Amazon River, near the Manaus City, Brazil, was computed, where the Negro River inflows into the Solimões River, from the North and the West, to form the Amazon River below this city. It is well known that in the Amazon's water system, the confluent tributaries usually have concentration difference in comparison with the mainstream, because of the humus in tropical rain forest, produced by tropic rains. The Negro River, however, is the largest black-water river in the world and also the largest left tributary of the river. The coarse yellow lines in this figure demonstrate the outline of computational domain. It is clear that the simulated depth-integrated contour is well coincident with the outline of black-water plume. In the computation, the original turbulence empirical constants of three depth-integrated models, suggested by their authors, are adopted and do not been changed never.   Figure 2 displays the digital map, on which the solved river reach with 56 short cross-river lines (i.e., NLrs=56) has been divided into 55 sub-reaches by the GUI of Q3drm1.0. It is notable that the cross-river lines within the riverside and island boundary have been redrawn in order to involve the configurations of two islands. Figure 3 and 4 illustrate the body-fitted non-orthogonal grids, generated by the Grid-Generator and plotted by the Grid-Browser of this software, with the resolution 227×18 on coarse grid and resolution 452×34 on the fine grid, respectively. The total length of the calculated river reach is up to 71.637km. In the generated two levels' grids, the nodal points on transversal grid lines are uniform. The tributary flows into the mainstream on the south riverside, with the numbers of nodal points from i=58 to 61 on the coarse grid. The two 'separated' islands start at (i=50, j=12) and (i=159, j=8), and ends at (i=75, j=12) and (i=195, j=8) on the same mesh, where the 'separated' means NICVS(2)>NICVE(1), i.e. 159>75. The flow direction is from the West and South to the East. The Grid-Generator generated two layers' grids, on which all of geometric data, necessary for the later calculation of flow and contaminant transport, should be stored and then can be read in by the Flow-Solver. It means that one CV on the coarse grid was divided into four CV on the fine grid. The bottom topography on fine grid, drawn by the Field Browser of Q3drm1.0, is presented by Figure 5. The variation of bottom topography will be considered during the calculation.

Flow and Side Discharge Solving
By using the developed Flow-Solver, the environmental flow and transport behaviors were simulated. In the Flow-Solver, Guass' divergence theorem, ILU (Incomplete Lower-Upper) decomposition, SIP (Strongly Implicit Procedure), PWIM (Pressure Weighting Interpolation Method), SIMPLE (Semi-Implicit Method for Pressure-Linked Equation) algorithm for FVA (Finite Volume Approach), under relaxation and multi-grid iterative method have been adopted. Firstly, the hydrodynamic fundamental governing equations were solved at the coarse grid and then at the fine grid, in the following sequence at each grid level: two momentum equations ( u -Eq. and v -Eq.), one pressure-correction equation ( ' p -Eq.), one concentration transport equation ( 1 C -Eq.), and two transport equations (i.e., the k -Eq. and  -Eq.; or k -Eq. and w -Eq.; or k -Eq. and  -Eq.), respectively.   Fig.6 to Fig.12. Fig.6 display the results, calculated by using   k model and drawn by the Field Browser, with a: flow pattern, b: streamlines, c: pressure field, d: concentration contours, e: k field and f:  field, respectively. Fig.6d illustrates that at the lower reach of the tributary outlet section, the pollutant plume well develops along the right riverside. The distribution graphics of the same physical variables and turbulent parameter k , calculated by turbulence  Fig.8b, ranges only from 1.936e-6 to 0.01162m 2 /s 3 ; however, the w value and  value range from 3.0524e-5 to 1.771/s 2 and from 0.547e-2 to 1.327/s in Fig.8c and Fig.8a respectively. Figs.9a, 9b and 9c display the 3D distribution graphics of effective viscosity eff  , while the depth-integrated turbulent eddy viscosity t  was calculated by using Eq.  Figs.11a and 11b show the comparisons of concentration profiles along the centers of CV at i from 1 to 456 and j=2 (i.e., along a curved line from the inlet to the outlet near the south riverside) and at i=350 and j from 1 to 38 (i.e., along a transversal section of i=350, which passes through the larger island in computational domain), calculated by three turbulence model closures, respectively. Fig.12a illustrates the comparisons between  , w and  along the curved line at j=3, and Fig.12b the comparisons of these three variables at i=350's transversal section. It is well known that the orders of magnitudes of  , w and  , adopted in the three turbulence models, have significant differences indeed.

The Development of Contaminant Plume at the Discharge Beginning
In order to well understand the development process of pollutant plume, the authors also performed a special calculation by using the closure of   k model for the case described as follows. Supposing the pollutant concentration through the tributary to equal zero firstly, and then, the concentration value instantaneously comes to 50mg/L at Time=0, while the flow-rates, either of the tributary or of the main stream, remain unchanged. Figs.13a-f display the development and variation of pollutant plume in the lower reach of tributary outlet section.

Discussions and Conclusions
It is well known that quasi 3D measurement data need to be obtained from depth-integrating real 3D measurement data, which is both rare and expensive indeed. Providing the closure function of multiple two-equation turbulence models not only enables users to make use of the new achievements of turbulence modeling theory, but also greatly elevates the credibility of the calculated results.
Because there is no an 'universal' closed model in turbulence modeling theory, providing several available turbulence closure models and letting users themselves to choose the suitable closure model(s) according to the solved problems probably become the best way to deal with the selection of turbulence closure models. The currently developed computing engines of advanced CFD software for the 'standard' 2D and 3D modeling in industry have widely adopted this advanced concept, such as the CFX, Phoenics and Star-CD software, developed by British companies, and the Fluent, Cfdesign and FIDAP software, developed by American companies, all can provide a number of two-equation turbulence closure model for users to choose and utilize.
Although the two-equations closure turbulence models are also very much an active research area, and new refined two-equation models with high-precision are still being developed, unfortunately, the existing quasi 3D simulation software and their corresponding mathematical models only can provide users with a depth-integrated two-equation turbulence   k model, appeared 40 years ago. Q3drm1.0 software, which can work on various Windows-based microcomputer, have overcome this hysteresis state and can provide multiple depth-integrated two-equations turbulence models for users. It is greatly convenient for users to select different turbulence models for closing the hydrodynamics fundamental governing equations and solving their problems in engineering.
Q3drm1.0 software focuses on the refined simulations of the steady and unsteady flow and temperature/contaminant transport problems in complicated computational domains with the strong ability to deal with different types of discharges: side-discharge, point-source discharge/point-sink, and area-source discharge from the slope along bank. However, only the study of side-discharge is presented in the paper.
At present, the turbulence k-ω model, just like the turbulence k-ε model, has become an industrial standard model and is widely adopted for most types of engineering problems. Therefore, the developed depth-integrated turbulence   k model and their numerical investigation and comparison with existing depth-integrated models, presented in this paper, are quite significant.
The distribution graphics of turbulent variable k are illustrated in Fig.7, which are calculated by three turbulence models, strongly vary in the computational domain, but still quite similar to one another. However, the distribution characteristics of  ,  and w , shown in Figs.8a, 8b and 8c respectively, are varying vary sharply and quite different from one another. In Figs.9a, 9b and 9c, the calculated effective viscosity eff  also varies strongly. Actually, the eddy viscosity changes from point to point in the computational domain, especially in the areas near the islands' boundaries and riversides. To solve the contaminant transport problems caused by side discharge, for example, the plume usually develops along a long and narrow region near riverside (see Fig.6d and Fig.13), where t  (or eff  ) varies much strongly indeed (see Fig.9). This means that t  cannot be roughly considered as an adjustable constant, and has to be calculated carefully and precisely by using suitable higher-order turbulence closure models with higher precision.
The concentration profiles along the south riverbank in Fig.11, either calculated by closures, or by w k closure, only have a quite small difference from one another. It means that three used depth-integrated turbulence models almost have the same ability to simulate plume distributions along riverbank. This conclusion also coincides with the previous research result of the first author that the depth-integrated two-equation turbulence models are suitable for modeling strong mixing turbulence [17]. The abilities of different depth-integrated two-equation turbulence models for rather weak mixing, also often encountered in engineering, should be investigated further.
It is notable that the order of magnitude of  is smaller than the order of magnitude of w , and much smaller than the order of magnitude of  . These three variables, however, all appear in the denominators of Eqs. (9), (14) and (2) for calculating turbulent eddy viscosity t  . For numerical modeling, the occurrence of numerical error is unavoidable, especially in the area near irregular boundary. A smaller numerical error, caused by solving  -Eq., for example, will bring on larger error for calculating eddy viscosity than the same error caused by solving other two equations (i.e., w -Eq. and  -Eq.). Doubtless, the elevation of the order of magnitude of the used second turbulent variable, which is reflecting the advance of turbulence closure modeling, provides a possibility for users to improve the computational precision. The insufficiency of traditional   k model may be avoided by adopting other turbulence models, appeared recently, such as the   k model and his revised version.
Two levels' grids in the current example, one coarse grid and one fine grid, were adopted. The simulations on these two level's grids can well satisfy the solved problem demand. If it is necessary, by changing the given number of grid levels at three in the software, modeling not only on coarse and fine grids, but also on finest grid can then be realized. The selected number of grid levels depends on modeler's requirements and the solved problems.
The solved depth-integrated concentration variable in present computation is the contaminant concentration difference between main stream and confluent tributary (50mg/L). In fact, other discharged contaminant indexes, such as BOD and COD, also can be considered as the solved variable. Q3drm1.0 software possesses the ability to solve simultaneously two concentration components in one calculation, produced by industrial and domestic discharges.
Q3drm1.0 software can be used to simulate, predict and analyze various scientific, engineering and application problems about river pollution, accidental discharge, water resources security, water environment protection, environmental engineering, environmental assessment and the comparison between different schemes for water supply and drainage and so on, which are closely related to flow, mixing and contaminant/waste heat transport.