Deng Yun1, Zhao Wenqian2, Li Jia3 and Luo Lin4
1Ph.D Candidate, 2,3,4 Professors
Sichuan University, the state key laboratory of high speed flow
Tel: 028-5407780; Fax: 5401558; E-mail: dengyun@fm365.com
Abstract: In this study the laterally averaged two-dimensional k–ε turbulent model is firstly tested in Ertan reservoir and the results are compared favorably with the observed measurements in site. Then it is applied to predict the thermal stratification for another reservoir. The forming, growth and erosion of the thermocline are simulated and the development of the thermal structure is analyzed all-around.
Keywords: thermal stratification, mathematical simulation, reservoir
The thermal stratification is one of important features, which has great impact on the ecological environment in the reservoir area and the downstream. How to simulate and predict the thermal structure for a huge reservoir has become a key problem for the reservoir design and management. The buoyancy due to the temperature difference greatly influences the flow field, which brings great difficulties to the computation. Jiang(2000) calibrated the turbulent viscosity in different direction for different variables with the observed data. Chu(1997) expressed the turbulent viscosity as the function of velocity gradient, thermal stratification parameter and turbulent length scale. The k–ε turbulent model with buoyancy modification has great advantages in stabilization and few parameters which is adopted by many researchers, such as Chen(1997) and Gerard(1988). For the simulation of a large reservoir in a long period, the cost of three-dimensional computation is too expensive. In the narrow reservoir that the variations of the variables in the transverse direction is not obvious, the laterally averaged two-dimensional model is a good choice which can well simulate the motion of the buoyant flow in the longitudinal and vertical directions and the development of the vertical thermal stratification. In this study, the laterally averaged two-dimensional k–ε turbulent model is verified with observed field data and applied to predict the thermal stratification in a huge reservoir.
2.1.1 The state equation
For the water in nature, the density can be expressed
(1)
The temperature difference is the buoyancy source. When the temperature difference is small, the buoyancy (gravity term) can be written as the linear function of the temperature difference
,
(2)
where
and
is the reference temperature and
density;
is the coefficient of thermal expansion.
The Boussinesq approximation is used in this study, which means that the buoyant effect is considered only in the gravity term.
2.1.2 The continuity equation
(3)
2.1.3 The momentum equations
(4)
(5)
2.1.4 The temperature equation
(6)
In the above equations u and w are the mean longitudinal
and vertical velocities respectively; T is the mean temperature; P
is the pressure; B is the river width;
, where
is the molecular viscosity and
is the kinematic eddy viscosity;
is the turbulent Prandtl number.
is specific heat of water;
is the solar radiation penetrating
into the water.
The standard k–ε turbulent model is used with modification by buoyancy forces. The
eddy viscosity
is computed from the equation
(7)
where k and ε are the turbulent kinetic energy and the kinetic
dissipation rate respectively.
The k and ε equations for the buoyant flow are:
(8)
(9)
in which
is the production of kinetic energy
from the mean flow and
is the production or destruction of turbulent kinetic energy by buoyant forces,
which are given by the expressions
,
(10)
The buoyancy tends to suppress turbulence and turbulent diffusion in the stably thermal stratification environment and boost up in reverse.
,
, the Prandtl number of the kinetic and kinetic dissipation rate, are given for the standard value of 1.0
and 1.3. The other constants are
:
The symmetry conditions for velocities are given on the free surface. The reservoir bottom and the dam face are taken as no-slip boundary. At the inflow and outflow boundary u is set as a constant and w as zero.
The bottom of reservoir and the surface of dam are taken as adiabatic. The inflow temperature is given a constant value and a zero gradient condition is set at the outflow boundary. The conditions on the free surface are much more complicated, since the heat exchange is the main heat source of reservoir. The heat fluxes through the water surface mainly include five parts
(11)
The net solar radiation
(12)
where
is incoming solar radiation;
(0.1) is the water reflection coefficient, and
(0.65) is the fraction of solar radiation absorbed at surface. The remaining
fraction of the solar radiation
is absorbed exponentially with
depth
(13)
where H is the water depth and
(0.5) is extinction coefficient.
The atmospheric longwave
radiation
(14)
where Ta is air temperature (Celsius degrees),
(=5.67×10-8W/m2·K4) is Stefan-Boltzmann constant,
and
is atmospheric emissivity
K (0.17) is the cloud height coefficient and C is the cloud cover ratio.
The back radiation from water to the air follows the same formulation with air temperature replaced by the water surface temperature Ts
(15)
and the emissivity
is fixed at 0.975
The evaporative loss (latent heat)
(16)
where es and ea
are the saturation vapor pressure above the water surface and the vapor pressure;
W is the air velocity 10 meters above
the water surface. The function of W is expressed
The conductive loss (sensible heat)
(17)
Zero gradient conditions for k andε are given on the free surface, the reservoir bottom, the outflow boundary and the dam surface. At the inflow boundary, k andε are computed by
,
(18)
where u0 and H0 are the water velocity and the water depth at the inflow boundary.
SIMPLE algorithm with staggered grid is used in this study (Patankar, 1980). Integration of governing equations over the control volume gives the discretized equations that are solved by TDMA (Tri-Diagonal Matrix Algorithm).
In order to verify the k–ε turbulent model for prediction of thermal stratification in the huge and deep reservoir, Ertan Reservoir is chosen as a test example. Then the model is applied to another huge reservoir, which is ongoing, for prediction its thermal structure.
Ertan reservoir locates in Southwest of China. The length is about 136 km, width about 400 m and maximum water depth about 190m with a huge storage capacity of 4.478×109 m3. It is a long, narrow, deep and huge reservoir, and is suitable to apply laterally averaged two-dimensional stratified model. The bottom of dam is at 1010m above the sea level. The maximum water level in the reservoir is 1200m in dry season. The intake of power station is located at 1128m with 10m high. In flood season, the superfluous water is discarded from the spillways and discharge orifices.
The temperature distribution in the reservoir is measured in January 14th of 2000, which is taken as initial conditions of temperature. Since the simultaneous velocity is not known, the velocity of non-buoyant flow is computed as the initial conditions. For lack of the simultaneous meteorological data, the monthly averaged meteorological data in 1999 at Shengli station is used and the long-time average water temperature at the end of reservoir upstream is taken as inflow boundary condition. The flow rate and water level are all monthly average shown in Fig. 1.
The reservoir area is discretized into 368×67 rectangle elements, the size varying from 10 to 400 meters in longitudinal direction and 2 to 3 meters in vertical direction. The fixed mesh is produced for each month with different water level.

Fig.1 Variation of discharge and water level with time
Numerical simulations from January to July of 2000 are performed and compared with the field observed data. Figure 2 shows the initial temperature distribution in January 14th. Figure 3(a)-(f) compares the computed and measured temperature distribution in the vertical cross section 700 meters away from the dam site upstream. It is surprised that they coincide fairly well especially in March and April under the conditions that the meteorological data and boundary conditions are not so accuracy. The results from May to July are not so good, due to the frequent reservoir operation (as shown in Fig. 1).

Fig.2
Temperature Distribution in January 14th
Figure 4 compares the computed outflow temperature and the temperature measured at the downstream hydrometric station. The result is similar with Fig. 3 that the lines of the computed and measured coincide very well from January to April and the computed temperature curve is much more smoothly than the measured. Before May, though the water level decreases fast, the flow rate varies smoothly and the water is all discharged from the intake of power station, so the outflow temperature rises steadily and is simulated well. The temperature downstream rises abruptly in late of May which conforms with the flood discharging from the tunnel spillway and the sequent frequent reservoir operation leads to fluctuating of the outflow temperature. The reservoir operation has great impact on outflow temperature, which can be used to control the outflow temperature.

Fig.3
Comparison of temperature vertical profile
¿
measured ——computed

Fig.4
Comparison of outflow temperature
In this model, only the solar radiation fraction and the extinction coefficient are needed to be determined. Through several trails, 0.65 and 0.5 are chosen for them. In fact they just appreciably influences the temperature of the water in 10 meters below the surface. Few parameters make the model very suitable for prediction.
The temperature model is applied to a huge reservoir, which is ongoing to construct. Because there are important fisheries downstream, the temperature of discharged water from reservoir is greatly concerned.
The reservoir is 200
km long and about 600 m wide. It is a typical river-like reservoir. The maximum
water depth is 220 m with a total storage capacity of 1.16×1010 m3.
The intake of its power station is located at 67 meters bellow the normal pool
level with 10 m high, from which water are almost discharged to the downstream
in normal year.
The reservoir area is divided into 527×71 rectangle elements, the size of which is the same as Ertan. It is observed from Ertan reservoir that the temperature tends to be uniform in February. So the inflow temperature in February is taken as the initial temperature in the whole reservoir area, and the velocity of non-buoyant flow as initial velocity. The corresponding meteorological data, inflow rate and water level are adopted. In order to eliminate the error of initial conditions, after the one year computation, the computation are repeated from the newly computed temperature and velocity in January until the temperature difference between the adjacent iterations is small enough, and then the temperature computed in the newest iteration is taken as the final result.
Figure 5 shows the
development of the temperature structure in four months, representing four
seasons, in the normal year. The temperature tends to be uniform in February
with the largest difference about 2℃ in vertical profile. In May with the increasing of
inflow temperature and heat flux from surface, the upper water got warmer fast
with the great temperature gradient forming at upper. But the water temperature
still keeps low near the bottom. The apparent thermal stratification forms. In
August the inflow and air temperature reach ultimate; there exists double
temperature gradient in the vertical profile. Due to the heat exchange between
water and air, a relatively small temperature gradient forms in a thin layer at
surface. Beneath it, the increasing flow rate enforces the water convection on
vertical section, so a warm and thermally uniform layer forms near the outlet.
The cold water in deep bottom is slowly warmed up and a thermocline layer exists
between the warm layer and the cold layer. In November, the air and the inflow
temperature decreases soon, the water loses heat to air and gets cold soon. With
cold-water gets heavier and goes down, the thermocline near surface disappears.
At the same time, the inflow gets cold and goes down along the bottom, which
makes the water near the bottom warmed up. The temperature difference declines
soon. With the weather gets colder, the temperature declines in whole vertical
profile until it tends to be uniform in the second February.
Figure 6 shows the
courses of air temperature, inflow temperature and the predicted water
temperature at surface and bottom at the dam site, and compared the outflow
temperature before and after reservoir built. The temperature at surface varies
soon, due to the direct heat exchange with the air, and the temperature
variation is small at bottom less than 4℃. Thermal stratification is apparent in the huge
reservoir. On the other hand, the temperature deferment of the huge reservoir is
obvious. February and October are the cross junctions. From March to September
after reservoir built, the outflow temperature declines 1.5℃ averagely, with the maximum of 3.5℃ in April. From November to the next January, the
temperature rises 2.8℃ than before with the maximum of 3.5℃ in December. The thermal storage
capacity also makes the maximum temperature difference in the whole year decline
from 12.4℃ to 9.9℃. Though the mean temperature in one
year is almost not changed, the temperature course changes much, which possibly
influences the ecological environment at downstream.

Fig.5
Temperature structure in different time
x—Distance from
upstream(m) y—water level(m)

Fig.6 Temperature variations with time(℃)
The mathematic model described in this paper could give realistic simulations of development, growth and erosion of the thermal stratification in the reservoir. Few parameters need to be calibrated. They are insensitive to the thermal structure for the deep reservoirs, which makes the model suitable for prediction.
The computed results in the two huge and deep reservoirs show that the thermal stratification exists almost in the whole year and displays the different temperature profiles in different seasons. The surface temperature is influenced directly by the climate conditions and varies fast. Reversely, the bottom temperature keeps cold in whole year and reaches the warmest in the autumn. The temperature deferment is obvious for the huge reservoir especially in the April and December. As to the outflow temperature, the temperature difference in one-year declines and its variation changes greatly, which may influence the ecological environment at downstream.
References
[1] Jiang Chun-bo, Zhang Qing-hai, Gao Zhong-xin. A 2-D unsteady flow model for predicting temperature and pollutant distribution in vertical cross section of a river. Journal of Hydraulic Engineering, Chinese. 2000, (9): 20-24.
[2] C.R. CHU and C.K. SOONG. Numerical simulation of wind-induced entrainment in a stably stratified water basin. J Hydraulic Research, Vol. 35, 1997, No.1.
[3] Gerard J. Farrell, Heinz G. Stefan. Mathematical modeling of plunging reservoir, J. Hydraulic Research. Vol. 26, No. 5, 1988.
[4] Chen Xiaohong, Liu Meinan, Lin Yanshan. Study on two-dimensional water quality distribution in reservoirs. Journal of Hydraulic Engineering, Chinese. 1997, (4): 9-16.
[5] PATANKAR, S.V.(1980). Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.
* Supported by the Key Program of the National Natural Science Foundation of China (Granted No. 59639240)