Effects of warm water inflows on the dispersion of pollutants in small reservoirs

ARTICLE IN PRESS

Journal of Environmental Management 81 (2006) 210–222 www.elsevier.com/locate/jenvman

Effects of warm water inflows on the dispersion of pollutants in small reservoirs Marı´ a C. Palancara,, Jose´ M. Arago´na, Fernando Sa´ncheza, Roberto Gilb a

Department of Chemical Engineering, Faculty of Chemistry, University Complutense of Madrid, 28040 Madrid Spain b Nuclear Safety Council, c/ Justo Dorado, 11, 28040 Madrid, Spain Received 27 July 2004; received in revised form 21 October 2005; accepted 24 October 2005 Available online 30 March 2006

Abstract The effects of the warm water discharged by a nuclear power plant (NPP) into a small reservoir are studied. A case study is presented (Jose´ Cabrera NPP—Zorita Hidra´ulica Reservoir) with experimental data of the reservoir stratification and predicted data of the dispersion of radioactive pollutants from operative or accidental releases. The vertical and longitudinal temperature profiles, electrical conductivity and transparency of the reservoir water were measured for an annual cycle. The results indicate that the continuous warm water discharge from the NPP causes permanent and artificial reservoir stratification. The stratification is significant within 1500 m upstream and 1000 m downstream from the warm water outfall. The pollutant dispersion has been predicted by using a flow model based on NT perfect-mixing compartments in series with feedback. The model parameter, NT, is calculated from the longitudinal diffusion coefficient. The prediction of pollutant dispersion by means of this model shows that the stratification slows down the vertical mixing in the whole water body, and reduces the reservoir volume that is effective for the dilution and dispersion of pollutants. This means that, in the case of a radioactive pollutant release, the reservoir radioactivity level could increase significantly. r 2006 Elsevier Ltd. All rights reserved. Keywords: Pollutant dispersion; Reservoirs; Warm water discharge; Radioactive pollutant

1. Introduction The thermal structure of a reservoir is the main regulator of almost all the physical, chemical and biological cycles (Wetang’ula, 2001). The variation of water density with temperature is the main cause for the thermal structure of reservoirs. The water density has a maximum value at a temperature within 0 and 3.94 1C. If the temperature is outside this range, the density decreases; this means that a vertical temperature gradient in a reservoir causes density gradients and buoyancy forces. These forces can impede the natural mixing processes in the mass of water due to wind and inflow effects. The buoyancy effects exist when the temperature of the upper layers of water is greater than the temperature of the deeper layers, i.e. the density is greater in the deeper layers. The stratification of the Corresponding author. Tel./fax: +34 91 394 4173.

E-mail addresses: [email protected] (M.C. Palancar), [email protected] (R. Gil). 0301-4797/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jenvman.2005.10.016

reservoir occurs when the buoyancy forces are greater than those of mixing. The stratification divides the mass of water into three characteristic layers: upper layer (epilimnion), intermediate layer (metalimnion) and deep layer (hypolimnion). The mass and heat transfer is practically impossible between these layers. In a stratified reservoir, the epilimnion allows abundant growth of algae, contributing to its eutrophication; hence the thermal stratification affects the quality of the water. In the case of wastewater releases, the stratification is a barrier for the mass and heat transfer and reduces the effective volume of the reservoir for the dilution and dispersion of pollutants with respect to the same reservoir in a homogeneous state. Temperature gradients and subsequent stratification occur in a reservoir by the heat exchange of water with the atmosphere and by natural or artificial inflows of water at a different temperature or salinity from the ones existing inside the reservoir. Only a small difference in temperatures is required to avoid the complete free circulation of the whole water column. The release of a warm stream from a

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

power plant can modify the water temperature (Holly et al., 2001), and produce stratification in the reservoir that receives the stream (Sinha et al., 1998). Both the magnitude of the intrusion and the formation of particular streamlines inside the reservoir depend on the volume and flow rate of the stream in relation to the reservoir volume and the main inflow rate (Bo-Ping et al., 2000). When there is a large enough discharge, a channel can be formed along the epilimnion that can remain along the whole reservoir (Ahlfeld et al., 2003). When the discharge occurs near the upper surface of the reservoir and the thickness of the inflow streamlines is smaller than that of the epilimnion, the intrusion stream can reach the dam discharge point. In this case, the dwelling time of the hypolimnetic water can be so long that the oxidative processes will increase the deoxyigenation of the hypolimnion. Whereas the natural stratification is a transitory and cyclic phenomenon, the stratification induced by the release of warm water from the nuclear power plant (NPP) refrigeration systems practically remains at a steady state, unless strong external forces destroy the stratification (e.g. a large cool water inflow). Ambient advection and diffusion processes control the mixing process in the far-field. Some models have been proposed to calculate the dispersion of routine and accidental releases into the surface waters in the far-field phase (Heling et al., 1999). The International Atomic Energy Agency (IAEA, 1985) clusters these models in three types: (1) numerical models that transform the basic governing equations of radionuclide dispersion into either finite differences or finite element form (Manson and Wallis, 2000), (2) box-type models that treat the entire body of water, or a section of the water body, as homogeneous compartments (Jacobsen et al., 1997) and (3) analytical models that solve the basic equations by describing the radionuclide transport (Fares, 2000; Gupta and Deshpande, 2004). The numerical models allow for relevant physical phenomena. These models are very complex and require knowledge of a great number of impoundment parameters. The box-type models compute average concentrations for each box (compartment) and set up transfer constants that serve to relate the variables to the adjacent compartments. Finally, the analytical models include major simplifications in order to make resolvable a complex analytical problem. The IAEA (IAEA, 1985) concludes that the simplest models, types 2 and 3, are usually satisfactory for NPP location purposes. Several box-type models can be applied to model impoundments. These models include completely mixed volumes, plug flow and partially mixed ponds, stratified reservoir, etc. The application of one or another of the models depends on the characteristics of both the reservoir and the power plant release. It is the objective of the present work to determine the influence of the inflow of the refrigeration cycle of a NPP on the dispersion and dilution of the release inside a small reservoir. Different methods are applied to calculate

211

vertical and horizontal dispersion coefficients in the reservoir mass of water. A conceptual model is proposed to calculate the pollutant dispersion in a permanently stratified reservoir. The dispersion in a non-stratified reservoir and in a permanently stratified reservoir is compared. A case study is applied using data of the warm water release from the Jose Cabrera NPP in the Zorita Hidra´ulica Reservoir. The paper is organized as follows: firstly, the hydrologic systems under study, the experimental methodology and the experimental results are described. Then, an experimental study of the water column structure (density, Brunt-Va¨isa¨la¨ frequency and vertical diffusion coefficient) is presented. A conceptual model, which serves to calculate the pollutant dispersion in a permanently stratified reservoir, is proposed and explained. The model is applied to calculate the pollutant distributions in two extreme cases, where the reservoir is in homogeneous state and where it is stratified as a consequence of the inflow. Finally, the predicted results are compared and discussed. 2. Reservoir description, experimental methodology and results 2.1. Reservoir description The Zorita Hidra´ulica Reservoir is used as the source of cooling water by the condenser of the Jose´ Cabrera NPP (nominal electrical power: 160 MW). The reservoir and NPP are located in the province of Guadalajara (Spain), as it is shown in Fig. 1. The reservoir is part of a system constructed on the Tagus River basin for different purposes: hydrological regulation, production of electrical energy and regional irrigation. A map of the Zorita Hidra´ulica Reservoir area is shown in Fig. 1. The reservoir has an extension of 7550 m and is limited by the Zorita and Bolarque Dams. The extraction is made by means of the reservoir outfall channel (downstream from a bridge, see Fig. 1) that drives the water to a hydraulic power plant (10.6 MWe). The top of this channel is the maximum level of water in the reservoir (600 m above Mediterranean Sea level). The channel has a depth of 6 m and is 12 m wide. The total volume of the reservoir is about 500 000 m3. There are two other reservoirs, the Almoguera and the Estremera downstream from the Zorita Hidra´ulica Dam. Both are used for irrigation, water supply for some villages and hydroelectric power plants. The Estremera Reservoir also supplies water to a fish hatchery. The water from the Zorita Hidra´ulica Reservoir flows towards the Tagus River bed and, for this reason, the reservoir has a rectilinear morphology for most of its extension. The reservoir makes a 901 bend only at a distance of 450 m upstream from the Zorita Dam (see Fig. 1). The cross-section of the reservoir bed is trapezoidal. Other characteristics of the reservoir are shown in Table 1. The Zorita Hidra´ulica Reservoir has no regulation function and therefore the inflow and outflow are

ARTICLE IN PRESS 212

M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

Fig. 1. Map of Spain with the six NPP currently in operation and details of the studied zone indicating the Tagus River, the Jose´ Cabrera NPP, dams, reservoirs (Zorita, Bolarque), and other details of interest.

Table 1 Characteristics of the Zorita Hidra´ulica Reservoir Maximum depth (m) Area of the surface layer (m2) Maximum capacity (m3) Useful capacity (m3) Maximum flow rate of flood (m3 s
1) Annual average inflow rate (m3 s
1) Annual average outflow rate (m3 s
1) Average width at 600 m above sea level (m) Average width at 592 m above sea level (m) Average dwelling time (h)

8 5.7  105 2.6  106 1.47  106 1800 10.35 10.24 92.3 10.4 39.5

practically the same. This means that the level of water in the reservoir is practically constant throughout the whole annual cycle. The values of water flow rate and level (datum) during an annual cycle are shown in Fig. 2. The NPP continuously extracts a flow of 6.67 m3 s
1 of cool water from the reservoir to carry out refrigerating operations. The cool water is pumped by means of two 1.37 m diameter inlet pipes whose centers are located at a depth of 3 m. The warm water produced is returned to the

reservoir by means of a surface channel that discharges downstream from the cool water inlet. The distances from the Bolarque Dam to the inlet and discharge points are 6500 and 6940 m, respectively. 2.2. Experimental methodology Twelve series of measurements were made monthly during the period November/97–October/98. The measured properties were as follows: temperature, electrical conductivity, water transparency and water layer thickness. These properties were measured at different depths and transversal cross-sections of the reservoir (Sa´nchez, 2000). The water temperature and electrical conductivity were measured by means of a conductivity meter with a probe, which includes a graphite cell and a Pt-100 RTD, connected to a 100 m electrical wire. The measures of electrical conductivity have an accuracy of 71 mS cm
1 and are given at the reference temperature of 20 1C. The accuracy of the Pt-100 is 70.1 1C. The measure of the water transparency was made through the Secchi depth with an accuracy of 70.1 m. The water layer thickness was measured by a sonar probe, with an accuracy of 70.1 m.

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

30 Datum

608

600 10

596 592

Conductivity (µS cm-1) 700 800

600

900

1000

0

-2 Depth (m)

604

20

500 Datum (m over sea level)

Flow rate (m3s-1)

Tagus R. inflow Tagus R. outflow

213

Temperature -4 Conductivity

0 0

50

100 150 200 250 Day of the annual cycle

300

-6

350

Fig. 2. Annual variation of the Tagus River inflow and outflow and datum of the Zorita Hidra´ulica Reservoir.

2.3. Experimental results Some representative examples of vertical profiles of temperature and electrical conductivity obtained in the runs are shown in Figs. 3 and 4. The measurements were made in a section of the reservoir located 607 m downstream from the discharge channel (near the bridge, see Fig. 1). The results show three main facts: (a) the temperature distributions prove that there are three layers in the water column, an upper and a lower layer that are almost uniform and an intermediate layer with a strong temperature gradient whose location depends on the date of measurement; (b) these layers are not detected by the vertical profiles of electrical conductivity; (c) it can be assumed from the data shown in Fig. 4 that there is a transversal homogeneity in the whole reservoir. The Froude number, calculated from the characteristics of the water column next to the Zorita Dam, is 36.57; this value is greater than 1/p. In addition, the dwelling time of the water in the reservoir is very short, no more than 160 h. These values of Froude number and dwelling time indicate that the reservoir cannot develop a natural stratification during its annual cycle. This has been experimentally confirmed by measuring several temperature profiles in the reservoir. As an example, the temperature profiles that are shown in Fig. 5 were obtained from a series of measures made in several sections of the Zorita Hidra´ulica Reservoir and within 4500 and 9000 m downstream from the

10

12

14

16 18 Temperature (°C)

20

22

24

Fig. 3. Vertical profiles of temperature and electrical conductivity in a sampling point located 607 m downstream from the NPP outfall channel. Date: 27/Feb/98.

0 Center Left side Rigth side

-2 Depth (m)

Governmental organizations and the companies that manage the reservoir and the NPP provided some additional information. The data provided are reservoir properties, such as evaporation loss, suspended solids, chlorophyll ‘‘a’’ concentration, and water temperature and flow rate discharged by the NPP. Additionally, the National Meteorological Institute provided information of the area, such as solar irradiation, cloudiness, air temperature, wind velocity and direction, air humidity and air pressure.

-8

-4

-6

-8 15

16

17

18 19 Temperature (°C)

20

21

Fig. 4. Vertical profiles of temperature at three points of a cross-section of the Zorita Hidra´ulica Reservoir located 607 m downstream from the NPP outfall channel. Date: 13/Oct/98.

Bolarque Dam. The distances in Fig. 5 refer to the first sampling point, which was located 4500 m downstream from the Bolarque Dam and 3500 m upstream from the NPP outfall channel. The last measured point was located at the Zorita Dam, which is located 9000 m downstream from the Bolarque Dam and 4500 m downstream from the first sampling point. The temperature in all the water column layers is constant from the Bolarque Dam to a distance of about 2000 m from the first sampling point. From 2000 m to the NPP outfall channel, located at 3500 m, the temperature of the upper layer (about 0–2 m deep) begins to increase and becomes constant downstream from the outfall channel. The temperature of the water

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

214

28 24 Temperature (°C)

Table 2 Energy inflows in the Zorita Hidra´ulica Reservoir within the surface area affected by the discharge

z=0m z=1m z=2m z=4m z=6m

20 16

Date

Solar irradiation (kJ m
2 d
1)

Inflow from NPP (kJ m
2 d
1)

8/May/98 16/Jul/98 13/Oct/98

21 000 29 000 16 000

40 000 78 000 51 000

12 inlet outfall 8 45 4

z=0m z=1m z=2m z=4m z=6m NPP inflow Tagus R. inflow

40 1000

2000 3000 Distance (m)

4000

5000

Fig. 5. Longitudinal profiles of temperature at various depths of the Zorita Hidra´ulica Reservoir. The abscise zero corresponds to a point located at 4500 m downstream from the Bolarque Dam and 3500 m upstream from the NPP outfall channel. Date: 13/Oct/98.

35 Temperature (°C)

0

30 25 20 15 10

winter

summer

5 0

50

100

150 200 250 Day of the annual cycle

300

350

Fig. 6. Annual variations of temperature of the NPP inflow, the Tagus River inflow and water at various depths of the Zorita Hidra´ulica Reservoir.

1000

Conductivity (µS cm-1)

layers at about 4–6 m depth remains constant in the whole reservoir. The behavior described above can be explained by considering that the temperature of the NPP liquid effluents is higher than the temperature of the Tagus River in the inflow area. This means an important heat input into the reservoir that can be higher than the contribution given by the solar irradiation. An example of typical values of both contributions is given in Table 2. The contribution of the NPP discharge can be up to four times the contribution of the solar irradiation. The annual variations of the NPP inflow temperature (inside the NPP outfall channel), the Tagus River inflow temperature (measured at 4500 m downstream from the Bolarque Dam) and the water temperature at various depths of the Zorita Hidra´ulica Reservoir (at a point located about 4 m in front the outlet of the NPP outfall channel) are shown in Fig. 6. During some periods of the year, as can be seen in Fig. 6, the water surface temperature is higher than the temperature of the surrounding air. Therefore, the energy losses of the water by sensible heat exchange can be the main contribution to the interchange with the atmosphere. Due to this fact, the surface layer of the reservoir follows the same temperature evolution as the NPP inflow warm water. In winter the heat exchange between the NPP discharge stream and the reservoir water becomes more important and the reservoir temperature, in all layers, is lower than the temperature of the NPP inflow warm water. The annual evolution of the temperature in the deepest water layers of the reservoir is similar to that of the Tagus River inflow. During summer, due to the low extinction coefficient value (Z ¼ 0:736 m
1 ) of the reservoir water and the low water layer thickness, the temperature of the deepest layers increases with time more significantly than the temperature of the Tagus River inflow.

z=0m z=1m z=2m z=4m z=6m

900

800

700

600 0

1000

2000 3000 Distance (m)

4000

5000

Fig. 7. Longitudinal profiles of electrical conductivity at various depths of the Zorita Hidra´ulica Reservoir. The abscise zero corresponds to a point located at 4500 m downstream from the Bolarque Dam and 3500 m upstream from the NPP outfall channel. Date: 13/Oct/98.

The longitudinal and vertical electrical conductivity profiles were measured at several reservoir points. It can be seen in the example shown in Fig. 7 that the electrical conductivity is practically constant throughout the reservoir. This fact indicates that the reservoir water

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

3. Determination of water column properties The experimental results obtained have served to study three important properties of the water column: water density, stability frequency (Brunt-Va¨isa¨la¨ frequency) and vertical diffusion coefficient. 3.1. Water density and stability frequency The water density profiles were calculated by means of Eqs. (1–4) (Chen and Millero, 1977), which are based on Kell’s equations (Kell, 1975), and serve to calculate the density and compressibility of the pure water, including the effects of solved salts (Millero et al., 1976; Chen and Millero, 1976). The accuracy of the equations is 72  10
6 g cm
3 within the following ranges of water properties: 0–0.06% of salinity, 0–30 1C of temperature and 0–180 bar of pressure. rP ¼

r0 1


P K

,

(1)

r0 ¼ b1 þ b2 T þ b3 T 2 þ    þ b7 T 6 þ ðb8 T þ b9 T 2 ÞS,

(2)

K ¼ c1 þ c2 T þ    þ c5 T 4 þ ðc6 þ c7 T þ c8 T 2 ÞP þ ðc9 þ c10 T þ c11 PÞS,

ð3Þ

S ¼ 0:9951437gT ,

(4)

P

0

where r and r are the water densities at current pressure and sea level pressure, respectively, K is a constant given by Eq. (3), P is pressure (bar), T is temperature (1C), and S is salinity, expressed by Eq. (4), where gT is the quantity of dissolved salt expressed as grams per kg of water. The stability of the water column was evaluated by means of the stability frequency (Brunt-Va¨isa¨la¨ frequency). This parameter expresses the strength of stratification in a fluid. It is used to evaluate the development of the thermocline in lakes, reservoirs and oceans (Jonas et al., 2003; Wunsch and Ferrari, 2004). The stability frequency, N 2, is defined by   qr g 2 , (5) N ¼
qz r where r is the water density, g the gravity acceleration and z the depth. Some authors (e.g. Wunsch and Ferrari, 2004) have used the potential density, s, to simplify the calculation of the stability frequency. In the case-study presented here it is not necessary to use s because the maximum reservoir depth is 8 m and the changes in density due to compressibility are negligible. Other authors (e.g. Garet and Ferron, 1996) independently evaluated the temperature and salinity influence. In our case, as it can be seen in

Figs. 3 and 7, the water reservoir and water discharge salinities are equal. Therefore, the stratification is only due to temperature differences. The experimental temperature profiles found in different sections of the Zorita Hidra´ulica Reservoir are well explained using two calculated characteristics: the water column density and the stability frequency. On the one hand, in the highest zone of the reservoir, from 1500 m upstream from the NPP discharge area to the Bolarque Dam, the profiles of temperature and electrical conductivity are practically flat. Consequently, the profiles of density and stability frequency are also flat. The stability frequency of the water column is low in all the measured points within this area; this is a typical characteristic of well-mixed systems. On the other hand, the density profiles in the area influenced by the warm NPP outfall show the high effects of the outfall temperature in the density of water. As an example, a density profile, calculated for a section near the Zorita Dam, is shown in Fig. 8. As the water temperature of the NPP discharge (about 15–25 1C) is higher than that of the reservoir, the water inflow has a smaller density than the reservoir water. The buoyancy effect of this warm water resists the vertical mixing action of wind. The continuous action of wind, combined with the action of the kinetic energy supplied by the discharge, develops a surface layer isotherm. The density of this layer is constant at a depth between 0 and 2 m in the example shown in Fig. 8, and is lower than the density of deeper water between 2 and 8 m in the example shown in Fig. 8. Therefore, this layer has the characteristic of an epilimnion. Below this layer, there is a transitional layer with a high hydrostatic stability due to a strong density gradient. Therefore, this layer is a metalimnion. The deepest layer presents the highest values of the water column density; it is the hypolimnion. There is a weak density gradient and the water density inside this

0.000

0.003

Stability (s-2) 0.006

0.009

0.012

0

-2 Depth (m)

composition is not affected by the use of water and the release of warm water from the NPP.

215

Density

-4

Stability -6

-8 997

998

999

1000

1001

1002

Density (kg m-3) Fig. 8. Vertical profiles of density and stability frequency in the column of water near the Zorita Dam. Date: 5/Feb/98.

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

layer is similar to the water density in the area near to the Bolarque Dam. In all the area between the Zorita Dam and the NPP discharge point, the structure of the water column is similar to the one described above and, consequently, the structure of this area is stratified artificially due to the NPP warm water discharge. The stability frequency curve in Fig. 8 clearly shows how the characteristics of the water column change abruptly at a depth of about 3 m due to the formation of the metalimnion. In fact, the stability frequency is a good index of the resistance to mass and heat transfer. Low values of stability frequency mean the mixing processes are strong, i.e. in the epilimnion and hypolimnion, while high values of the stability frequency are indicative of poor mixing between different layers of the water column, i.e. in the metalimnion. 3.2. Vertical diffusion coefficient The calculation of the vertical diffusion coefficient depends on the structure of the water column. For the water layer between the upper surface and the depth where the influence of wind and the surface heat input are negligible, the expression relating this coefficient with depth is   ra C D 3=2 u3a e z ¼ gm , (6) r N 2 kh where ez is the vertical diffusion coefficient (m2 s
1), gm is the fraction of the wind kinetic energy that is added to the buoyancy flow at each depth, ra is air density, CD is friction coefficient, ua is wind speed at 10 m height above the water surface, N2 is the stability frequency, k is the von Ka´rma´n constant (set at 0.41), and h is depth. Kullemberg (1973) has proved that Eq. (6) is adequate to explain the vertical diffusion of a tracer in a lake. The coefficient CD depends on a number of factors (Amorocho and de Vries, 1980; Hicks, 1972; Hicks et al, 1974; Wu, 1969). Nevertheless, a mean value of 1.3  10
3 may be used for engineering purposes, as has been recommended by Fischer (Fischer et al., 1979). The value of the mixing efficiency, gm, was set at 0.12 (Ivey and Imberger, 1991). Assuming that the whole reservoir is in equilibrium with the external energy inputs (Fischer et al., 1979), the vertical diffusion coefficient for deep water can be calculated by ez ¼

CH 2 , T m Sc

Tm ¼

Ep dE p Pe As þ dt

(7)

stability and Imberger and Patterson (1981) have set C at a constant value of 0.048, H is the reservoir depth, Tm is the mixing time, calculated by Eq. (8), Ep is the water potential energy, Pe is the wind power per surface unit, As is the surface of the reservoir, Sc is the stability and Dr is the difference in density between the surface and bottom layers. An example of the warm water influence on the water column structure is shown in Fig. 9, where the variation of the vertical diffusion coefficient with depth in two sections of the reservoir is shown. One section is located downstream, near the Zorita Dam; a zone clearly influenced by the warm NPP outfall. The other section is at 1700 m upstream from the NPP outfall channel; a zone not influenced by the warm NPP outfall. In all the measured points, the highest value of the vertical diffusion coefficient was found in the upper water layer (see, for example, the upper part of the two curves in Fig. 9). This means that the surface layer receives an important energy contribution from wind that directly produces mixing. The eddy energy induced by wind is gradually dissipated with depth because it acts against the gradient of density. The result is that the vertical diffusion coefficient decreases with depth. In the area not influenced by the warm water release, the decrease of this coefficient is gradual to the bottom of the reservoir (e.g. the curve on the right of Fig. 9). In the area where the warm water influence is significant, the diffusion coefficient abruptly decreases at the depth corresponding to the thermocline. At this depth, the stability frequency is very high (in fact, its value is a maximum) and the dissipation of the available kinetic energy is important; the main consequence is that the mixing process is strongly impeded. In the upper layer water, the vertical diffusion coefficient ranges between 10
4 and 10
5 m2 s
1. The upper limit corresponds to the points located in the zone not influenced by the NPP discharge; the lower limit corresponds to the points situated in the zone influenced by the discharge. The 0 -1 -2 Depth (m)

216

-3 -4 -5 -6

,

zone influence zone not influence

(8) -7 -8

dr H , Sc ¼ dz Dr

(9)

where ez is the vertical diffusion coefficient expressed in m2 s
1, C depends on the reservoir shape, stratification and

10-8.0

10-7.0

10-6.0 10-5.0 Diffusion coefficient (m2 s-1)

10-4.0

10-3.0

Fig. 9. Profiles of vertical diffusion coefficient next to the Zorita Dam and at 1700 m upstream from the NPP outfall channel. Date: 5/Feb/98.

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

difference between both values is due to the different water densities in both areas (999.505 kg m
3 in the noninfluenced area and 998.292 kg m
3 in the influenced area). In the influenced area, the low water density of the upper layer causes buoyancy forces that act against the up-todown mixing. These forces are lower in the not influenced area, in which the stability frequency is lower and the vertical diffusion coefficient is higher. Although the annual variation of the diffusion coefficient is affected by solar irradiation, the effect of the NPP warm water temperature on this coefficient is stronger. To illustrate this, the annual global irradiation of sun in the reservoir area and the NPP inflow temperature are shown in Fig. 10a. The vertical diffusion coefficient, calculated for the upper layer, thermocline and bottom layer, and the thermocline depth, at a point located about 4 m in front the outlet of the NPP outfall channel, are shown in Fig. 10b. Comparing the lines of Fig. 10b with the evolution of the NPP inflow temperature, Fig. 10a, it can be seen that, as the temperature of the NPP inflow increases, the value of the dispersion coefficient decreases. The effect is the same for every depth, the NPP inflow temperature being the most important effect on the value of the coefficient in the thermocline layer. It also becomes appreciable that when the thermal contribution of the discharge is the highest, around day 200, the thermocline depth is lower. This is a consequence of the highest heat transfer and the mixing action of wind. This influence is higher than the influence of solar irradiation, which is the highest around day 150,

40000

26

22 20

20000 18 16

10000

14 (a)

10-1

Diffusion coefficient (m2 s-1)

NPP inflow temperature (°C)

24

NPP inflow 30000

Upper layer Botton layer Thermocline

10-2

Thermocline depth 2.0

10-3 10-4

4.0

10-5

6.0

Depth (m)

Global irradiation (kJ m-2 d-1)

Irradiation

10-6 8.0

10-7 (b)

10-8 0

50

100 150 200 250 Day of the annual cycle

300

10.0

217

Fig. 10a. Three examples of the energy inflows into the Zorita Hidra´ulica Reservoir, in spring, summer and fall, are shown in Table 2. 4. Modeling the dispersion and dilution of a pollutant in the reservoir The model, which is described in detail in the following paragraphs, is useful to calculate the longitudinal dispersion from experimental measurements and to predict the concentration–time distributions upstream and downstream from the NPP discharge point. 4.1. Modeling overview The analysis of the water column structure shows that the reservoir is subjected to artificial stratification due to the NPP warm water discharge. The temperature range of the release is about 15–25 1C. The diffusion coefficient at the thermocline has values ranging between 10
6 and 10
7 m2 s
1; which are values close to those of the molecular diffusion. Therefore, the thermocline acts as a barrier for both the mass and energy transfer. Two direct consequences are that a pollutant release travels into the epilimnion of the water column and the actual reservoir volume for the dispersion and dilution of potential pollutant releases does not correspond to the total reservoir volume but to the volume of the epilimnion. Moreover, the stratification occurs upstream as well as downstream from the discharge; this implies that the water discharged from the NPP flows into both sides of the reservoir, upstream and downstream from the discharge point. This means that a discharge could not only affect the areas located downstream from the discharge but also the areas located upstream. To quantify this fact, the model represented in Fig. 11 is proposed. The far field is treated here by assuming that several perfect-mixing compartments of equal volume connected in series compose the epilimnion of the water column. This model has been used to describe the flow through real vessels (Levenspiel, 1979) and serves here to describe the behavior of the hydrological system. The departure hypothesis is that the concentration of a solute in a given point of the reservoir can be calculated by the equation that describes the output of a system composed of NT perfect-mixing compartments. The model assumes that all the compartments have the same volume and the total volume of the model is equal to the total volume of the real system. The supposition of ‘‘perfect-mixing’’ implies that the concentration and temperature inside each compartment are uniform and equal to those in the compartment output.

350

Fig. 10. (a) Annual variation of the global irradiation and the temperature of the NPP outfall. (b) Annual variations of the vertical diffusion coefficient, and thermocline depth.

4.2. Description of mathematical model The model proposed to represent the epilimnion is based on a series of perfect mixing compartments of equal

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

218

QD CD Inflow from the NPP QR CR

QR-1 CR-1

QB CB0 1 Inflow from Bolarque

Qn-2 Cn-2

QS+1 CS+1

QS CS

n-1

n

n+1

QS+m CS+m n+m Output to Tagus

Qn-1 Cn-1

Upstream area influenced by the inflow from the NPP

Discharge area

Downstream area

Fig. 11. Model of perfect-mixing compartments in series.

volume. The model (see Fig. 11 for details and symbols) is applied under the following assumptions: (1) the water from the NPP flows into compartment n with a flow rate QD and a concentration CD; (2) compartment n has a second inflow of water that comes from compartment n
1, located immediately upstream from it, with a flow rate Qn
1 and a concentration Cn
1; (3) compartment n has two outflows, one is a feedback QR to the upstream compartment n
1 and the other one is a stream with flow rate QS to the downstream compartment n+1; and (4) there is no downstream water feedback from compartment n+1 to compartment n; (5) the reservoir area located upstream from the NPP discharge point can be considered as n
1 compartments in series (compartments 1 to n
1 in Fig. 11). Compartment 1 has two inputs; one is the water inflow coming from the Bolarque Reservoir, QB, the other corresponds to the feedback from compartment n
1 originated by the inflow of the discharge, QR
1. Compartment 1 has only one output; there is no feedback to the Bolarque Reservoir because the Bolarque Dam is an impassable barrier. Compartments 2 to n
1 have two inputs and two outputs each. The inputs correspond to the normal stream from the preceding compartment and the feedback from the subsequent compartment. The outputs correspond to the normal water flow to the subsequent compartment and the feedback to the preceding compartment, (6) the reservoir area downstream from the NPP discharge point can be considered as m compartments in series (n+1 to n+m in Fig. 11) with one input and one output each because, in this zone, there is no feedback. For a compartment n, the mass conservation balance of a solute in unsteady state is d ðV n C n Þ þ V n lC n , (10) dt where Cn
1 and Cn are the input and output pollutant concentrations in compartment n, respectively, Qe and Qs are the input and output flow rates of compartment n, respectively, Vn is compartment volume and l is the kinetic constant of a first-order non-conservative process. Eq. (10) is also valid for radioactive pollutants, in which case the concentration would be substituted Qe C n
1 ¼ Qs C n þ

by the specific activity and the kinetic constant by the radioactivity decay constant. The general balance, Eq. (10), for a non-conservative pollutant (e.g. a radioactive compound) in compartment n is valid for the section of the epilimnion located downstream from the discharge, i.e. in the zone where there is no feedback. Compartments n
1, n
2, etc. describe the section of the epilimnion located upstream from the discharge, i.e. in the zone with feedback. The mass balances for the pollutant in these compartments are different from the balance in compartment n described above. For example, for compartments n and n
1, the balances are given by QD C D þ Qn
1 C n
1 ¼ QR C R þ QS C S d þ ðV n C n Þ þ V n lC n , dt

ð11Þ

Qn
2 C n
2 þ QR C R ¼ QR
1 C R
1 þ Qn
1 C n
1 d þ ðV n
1 C n
1 Þ þ V n
1 lC n
1 . ð12Þ dt The model proposed predicts the concentration–time distributions of a pollutant at different distances from the NPP discharge point. The model input data are: Pollutant release characteristics (radioactive compounds, shape of the release: pulse, step, etc.), flow rate or volume and concentration of the pollutant, hydraulic and hydrological characteristics of the reservoir (change of capacity, surface section and width with depth, and variation of depth with distance from the beginning of the influence zone, thermocline depth and actual water flow rate). The parameter estimated is the number of compartments within each area influenced by the NPP release (i.e. n compartments upstream and m compartments downstream from the inflow point, see Fig. 11). 4.3. Determination of number of compartments, N T The total number of compartments in the epilimnion, NT, must be known to solve numerically Eqs. (10–12). The number of compartments, NT, has been calculated by Eq. (13), which has been obtained by combining

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

the dispersion and the tanks in series models proposed by Levenspiel (1979) NT ¼

u2m tm , 2DL

(13)

where um and tm are the mean velocity and residence time of the water, respectively, and DL is the longitudinal dispersion coefficient. The calculation of NT by Eq. (13) guarantees that the solution of the model with NT compartments reproduces the solution of the real dispersion equation describing the transport. The experimental determination of DL for a reservoir can be made by means of either dye or radioactive tracers. However, in the present case, the presence of the NPP made this type of determination impossible, given the fact that the sensitivity of people living around this type of facility is very high. Then, the use of colored tracers in a populated area would create social alarm (although the tracers commonly used are innocuous for the biota) and the use of radioactive isotopes—‘‘invisible’’ tracers—would increase the water activity in the areas located downstream from the reservoir. So, the alternative followed to calculate DL was to estimate it from the heat balance and following the hypothesis of Okubo (1980), which assumes that the horizontal turbulent diffusion coefficient depends on the scale of the process but not on the properties of the considered system (temperature, concentration of solved or suspended solids, viscosity, etc.). To estimate DL a two-dimensional model is used with the assumptions of constant volume, transversal homogeneity and horizontal isotherms. The reservoir is divided into cells of uniform size. The general expression of the energy conservation is jAccumulationj

Internal Absorption ¼ jDiffusionj þ of Radiation Horizontal
jConvectionj þ . Advection

ð14Þ

When the cell is within the surface layer, a term of net heat flow into the air–water interface is added to Eq. (14). If the cell is within the zone influenced by the NPP discharge inflow, a term of heat input due to the discharge is added to Eq. (14). The energy terms considered in the air–water interface are: net solar irradiation, net long wave radiation, latent heat of vaporization and sensible heat flux. The light absorption in the surface layer and in the whole water volume is given by Eqs. (15) and (16), respectively: R0 ¼ ð1
rf ÞbRsol ,

(15)

Rz ¼ ð1
rf Þð1
bÞRsol e
Zz ,

(16)

where Rsol is incident solar irradiation, rf is albedo, b is the factor of absorption surface, which was set at 0.4 (Dake and Harleman, 1969), and Z is the total extinction coefficient, which was calculated by the correlation of

219

Osgood (1984) and the experimental data of the Secchi depth. The net long wave radiation, back radiation, is defined as the difference between the long wave radiation emitted by the water, calculated by the equation of Stefan–Boltzmann, and the long wave radiation emitted by the atmosphere, given by Ra ¼ ð1
ra Þss ea T 4a ,

(17)

where ra is reflection of infrared radiation in the water surface, which was set at 0.03, ss is the Stefan–Boltzmann constant, Ta is the absolute air temperature and ea is the atmospheric emissivity, which was calculated by the correlation proposed by Idso and Jackson (1969). A large number of formulations have been proposed to calculate the latent heat of vaporization. The best fit was obtained here by the formulation proposed by Neuwirth (1974) H e ¼ f 0 ðua Þðew
ea Þ,

(18)

0

where f ðuÞ ¼ 3:72 þ 3:06ua ; ua is the wind speed at a meters above the water surface and ew and ea are vapor pressure on the water surface and in air at ‘‘a’’ meters from the water surface, respectively. The interchange of sensible heat in the surface of the reservoir is given by H c ¼ ra C p f ðuÞðT s
T a Þ,

(19)

where ra, Cp and Ta are density, specific heat coefficient for isobaric conditions and air temperature, respectively, f ðuÞ is a function of wind speed, f ðuÞ ¼ cua (Hondzo and Stefan, 1993). The coefficient c depends on the reservoir surface; according to Ford and Stefan (1980) it was calculated by c ¼ 24 þ lnðAk Þ, where Ak is the reservoir surface (km2). The mixing between the reservoir water and the outflow from the NPP in the discharge element (compartment n, see Fig. 11) is treated in the energy conservation equation, Eq. (14), by adding a term of heat input due to the discharge and by assuming that the inflow causes a reservoir water entrainment at the inflow depth and a mixing of the entrained water with the inflow. To estimate the temperature in the compartment n, it is necessary to know the amount of water entrained. The calculation method of amount of water entrained is poorly known in the modeling of reservoirs. Therefore, a relative value of the entrainment with respect to the inflow should be assumed. It was assumed here a flow rate ratio of 0.5, as recommended in RSTEMP (1972) for inflows at the water surface, which is the case of warm water discharge studied here. This value means that the entrainment is 50% of the inflow and, consequently, the total effective inflow is 1.5 times the actual stream inflow from the NPP. Finally, it is assumed that the discharge compartment is perfectly mixed and all the variables are known (by measure or PNN information) with the exception of the feedback flow (QR in Fig. 11) and the input of cool water to this compartment (Qn
1 in Fig. 11). These two flows are

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

220

estimated by applying the mass and energy conservation equations to the discharge compartment.

coefficient as a function of geometrical parameters DL ¼ 0; 011

4.4. Estimation of dispersion coefficient Some mean values of the longitudinal dispersion coefficient calculated for the Zorita Hidra´ulica Reservoir are shown in Fig. 12. The value of DL was estimated for an area located around a reference point 3500 m upstream from the NPP outfall channel. The value of DL is almost constant in all the studied areas except in the zone influenced by the warm discharge, at 3500 m in Fig. 12. The artificial epilimnion is formed between 3000 and 4000 m (see Fig. 12). In this area of influence, the dispersion coefficient decreases, upstream and downstream from the discharge point. The values of this coefficient range between 8 and 4 m2 s
1. The first value of DL is only for the point of discharge and is a large value that can be explained by the kinetic energy contribution introduced by the discharge. The reservoir area that is not influenced by the warm water released from the NPP presents similar hydrologic characteristics as the area where the Tagus River enters the Entrepen˜as Reservoir, another reservoir that was studied in a previous work (Sa´nchez, 2000). It was possible to make tracer runs in this zone of the Entrepen˜as Reservoir and to calculate the dispersion coefficient from the tracer response curves obtained. As this zone has similar dispersive and diffusive characteristics as the non-influenced zone of the Zorita Hidra´ulica Reservoir, the longitudinal dispersion coefficient obtained by heat balance (Okubo, 1980) was compared with the one determined from tracer data in the Entrepen˜as Reservoir. The differences in geometric parameters were corrected by the correlation proposed by Fischer et al. (1979) that serves to calculate the dispersion

Dispersion coefficient (m2 s-1)

10

8

6

4

2

0 0

1000

2000 3000 Distance (m)

4000

5000

Fig. 12. Mean value of the dispersion coefficient in the Zorita Hidra´ulica Reservoir. The abscise zero corresponds to a point located at 4500 m downstream from the Bolarque Dam and 3500 m upstream from the NPP outfall channel.

u2m W 2 , u d

(20)

where um is the mean velocity, W is the bed width, u* is the friction velocity, given by u ¼ ðg d Se Þ1=2 , where d is the mean depth, g is the acceleration of gravity and Se is the water surface slope. Eq. (20) was used to allow for the differences in geometrical characteristics between both reservoirs. The DL of the Entrepen˜as Reservoir was then corrected for the geometrical characteristics of the Zorita Hidra´ulica Reservoir. The corrected final value of DL was 1.6 m2 s
1, which is of the same order of magnitude as was obtained from Eqs. (14–20) for the non-influenced area of the Zorita Hidra´ulica Reservoir, i.e. about 4 m2 s
1 (see Fig. 12). 5. Discussion on simulation results The model equations, Eqs. (10–12), were implemented in a computer code that uses the method of finite time increments. The code includes the number of compartments, geometrical expressions to calculate the volume of each compartment and hydraulic and hydrologic characteristics. The code output gives the evolution of the pollutant concentration in the area of discharge and near the Zorita Dam. A numerical simulation was made by assuming a discharge in the form of an instantaneous impulse and the parameters that are shown in Table 3. For these conditions, the number of perfect-mixing compartments considered in the model is 64 and the volume of each one is 8214 m3. The time step was set at 300 s in the numerical integration of Eqs. (10–12); this is a time smaller than the residence time of each compartment. The results obtained are shown in Figs. 13 and 14. It can be clearly seen from these results that the artificial stratification affects the specific activity. It is higher than the activity calculated with the assumption of no stratification (for example, it is 42% higher at time zero, Fig. 13). This fact is not important if the specific activity is low enough, but it can be very important when the specific activity of the discharge is relatively high. In this case, the artificial stratification affects not only the pollutant concentration but also the time necessary for the response to pass through a given section. The nuclear safety regulations require the adoption of measures and decisions to minimize the impact of the discharges, for example the temporary prohibition of water consumption from the reservoir, or forbidding some recreational activities while the specific activity is higher Table 3 Conditions used in the simulation Radionuclide Activity release (Bq) Tagus River flow rate (m3 s
1) Datum (m over sea level)

H-3 3.7  1010 44.23 599

ARTICLE IN PRESS M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

Specific activity (Bq m-3)

5E+06 With stratification Without stratification

4E+06

3E+06

2E+06

1E+06

0E+00 0.0

0.1

0.2

0.3

0.4

0.5 0.6 Time (h)

0.7

0.8

0.9

1.0

Fig. 13. Distribution of pollutant concentration in the proximity of the NPP outfall channel.

Specific activity (Bq m-3)

3E+05 With stratification Without stratification 2E+05

1E+05

0E+00 2

3

4

5

6 7 Time (h)

8

9

10

11

Fig. 14. Distribution of pollutant concentration at a point near the Zorita Dam.

than a certain value. Assuming, for example, that the limit of specific activity permissible in the reservoir is 1.5  105 Bq m
3, the curves of Fig. 14 show that, neglecting the stratification, the prediction model gives values of the specific activity near the Zorita Dam that apparently do not reach that limit; consequently, no special safety measures are necessary. However, the real value of the specific activity is the one predicted by the model considering that the stratification exists; in such a case, from the lines of Fig. 14, it can be seen that the specific activity is higher than the given limit and some safety measures must be taken. 6. Conclusions The influence of the continuous warm water discharge from a NPP in the dynamic behavior of a small reservoir was studied. Several measures of water temperature, electrical conductivity and transparency were made at

221

different cross-sections of the Zorita Hidra´ulica Reservoir (Spain) throughout an annual cycle. The profiles of water density, stability frequency (Brunt-Va¨isa¨la¨ frequency) and vertical diffusion coefficient were calculated from the experimental measures as well as from other complementary information (meteorological, hydrological and geometrical data). These profiles were used to characterize the water column structure. The analysis of the results indicates that the warm water discharge causes an artificial stratification that exists during the whole time of the operation of the NPP and is not destroyed by normal meteorological or flow regime factors. This stratification produces a decrease in the vertical diffusion coefficient in the thermocline. This means that this coefficient can reach values close to the molecular diffusion and that the vertical mixing in the water column is impeded. The effects of the artificial stratification are significant from the outfall point up to 1500 m upstream from the outfall and to the reservoir dam, located 1000 m downstream from the outfall. When the longitudinal dispersion coefficient cannot be calculated by other methods, this parameter can be calculated by means of the general expression of the energy conservation equation. The values of this coefficient obtained by this method are close to those obtained by others methods. A box-type model, based on a series of uniform mixing compartments, has been successfully defined and applied. The model is able to include the feedback of the discharge that causes a back-mixing upstream from the outfall and to predict the artificial stratification. The prediction of stratification is important for decision makers in water management and nuclear safety areas since the specific activity and the time needed to pass a radioactive pollutant release is significantly greater in a stratified reservoir than in the same reservoir with a non-stratified structure.

Acknowledgments The study has been carried out with financial support from the Spanish Nuclear Safety Council (CSN).

References Ahlfeld, D., Joaquin, A., Tobiason, J., Mas, D., 2003. Case study: impact of reservoir stratification on interflow travel time. Journal of Hydraulic Engineering 129 (12), 966–975. Amorocho, J., de Vries, J.J., 1980. A new evaluation of the wind stress coefficient over water surfaces. Journal of Geophysical Research 85, 433–442. Bo-Ping, H., Argenmel, J., Garcı´ a, J.C., Cormena, M., Roura, M., Dolz, J., Straskraba, M., 2000. The thermal structure of Sau Reservoir (NE: Spain): a simulation approach. Ecological Modelling 125, 109–122. Chen, C.T., Millero, F.J., 1976. Specific volume of sea water at highpressures. Deep-Sea Research 23, 595–612. Chen, C.T., Millero, F.J., 1977. The use and misuse of pure water PVT properties for lakes waters. Nature 266, 70–71.

ARTICLE IN PRESS 222

M.C. Palancar et al. / Journal of Environmental Management 81 (2006) 210–222

Dake, J.M.K., Harleman, D.R.F., 1969. Thermal stratification in lakes: analytical and laboratory studies. Water Resources Research 5 (2), 484–495. Fares, Y.R., 2000. Dispersion processes in reservoir and natural stream. Proceedings of the Institution of Engineers 214 (6), 857–866. Fischer, H.B., List, E.J., Koh, R.C.Y., Imberger, J., Brooks, N.H., 1979. Mixing in Inland and Coastal Waters, first ed. Academic Press Inc., San Diego (Chapter 6). Ford, D.E., Stefan, H.G., 1980. Thermal predictions using integral energy models. Journal of Hydraulic Division, ASCE 106 (HY1), 39–55. Garet, A.E., Ferron, B., 1996. The effects of differential vertical diffusion of T and S in a box model of the thermocline circulation. Journal of Marine Research 54, 847–866. Gupta, S.K., Deshpande, R.D., 2004. An insight into the dynamics of Lake Nainital (Kumaun Himalaya, India) using stable isotope data. Hydrological Journal 49 (6), 1099–1113. Heling, R., Raskob, W., Popov, A., Zheleznyak, M., 1999. Overview of hydrological dispersion module-HDM of RODOS (1999) RODOSWG4-TN (99) 18. Report from the Comission of the European Comunities. Hicks, B.B., 1972. Some evaluations of drag and bulk transfer coefficients over water bodies of different sizes. Boundary–Layer Meteorology 3, 201–213. Hicks, B.B., Drnklow, R.L., Grauze, G., 1974. Drag and bulk transfer coefficients associated with a shallow water surface. Boundary–Layer Meteorology 6, 287–297. Holly Jr., F.M., Bradley, A., Rocco, B., 2001. Thermo-hydrodynamic modeling in management of nuclear energy production. Proceedings of IAHR XXIX Congress, Theme B: Environmental Hydraulics and Ecohydraulics, Beijing, China, September 16–21, pp. 409–416. Hondzo, M., Stefan, H.G., 1993. Lake water temperature simulation model. Journal of Hydraulic Engineering ASCE 119 (11), 1251–1273. Idso, S.B., Jackson, R.D., 1969. Thermal radiation from the atmosphere. Journal of Geophysical Research 74 (23), 5397–5403. Imberger, J., Patterson, J.C., 1981. A dynamic reservoir simulation model DYSRESM—5. In: Fisher, H.B. (Ed.), Transport Models for Inland and Coastal Waters. Academic Press Inc., New York, pp. 310–361. International Atomic Energy Agency (IAEA), 1985. Hydrological dispersion of radioactive material in relation to nuclear power plant sitting. Safety Series no 50-SG-S6. Ivey, C.N., Imberger, J., 1991. On the nature of turbulence in a stratified fluid, Part I: the energetics of mixing. Journal of Physical Oceanography 21, 650–659. Jacobsen, J.L., Madsen, H., Harremoe¨s, P.A., 1997. A stochastic model for two-station hydraulics exhibiting transient impact. Water Science Technology 36 (5), 19–26.

Jonas, T., Stips, A., Eugster, W., Wu¨est, A., 2003. Observations of a quasi shear-free lacustrine convective boundary layer: stratification and its implications on turbulence. Journal of Geophysical Research 108 (C10), 26.1–26.15. Kell, G.S., 1975. Density, thermal expansivity and compressibility of liquid water from 0 1C to 150 1C: correlations and tables for atmospheric pressure and saturation reviewed and expressed on 1968 temperature scale. Journal of Chemical and Engineering Data 20, 97–105. Kullemberg, G., 1973. An experimental study of diffusion characteristics in the thermocline and hypolimnion regions of Lake Ontario. Proceedings of the 16th Conference on Great Lakes Research, pp. 774–790. Levenspiel, O., 1979. The Chemical Reactor Omnibook, first ed. OSU Book Stores, Corvallis, OR (Chapter 66). Manson, J.R., Wallis, S.G., 2000. A conservative, semi-Lagrangian fate and transport model for fluvial systems—I. Theoretical development. Water Research 34 (15), 3769–3777. Millero, F.J., Gonza´lez, A., Ward, G.K., 1976. Density of seawater solutions at one atmosphere as a function of temperature and salinity. Journal of Marine Research 34, 61–93. Neuwirth, F., 1974. Grundlagen fu¨r die beurteilung der wa¨rmebelastung. Aufl. 1. Okubo, A., 1980. Diffusion and ecological problems. In: Krickeberg, K., Levin, S.A. (Eds.), Biomathematics. Springer, Berlin, pp. 11–15. Osgood, R.A., 1984. A 1984 study of the water quality of 43 metropolitan area lakes. Publication no 10-84-172, Metropolitan Council, St. Paul, MN. RSTEMP. 1972. Reservoir temperature stratification. User’ Manual, Hydrologic Center US Army Corps of Engineering, USA. Sa´nchez, F., 2000. Transporte y dispersio´n de contaminantes en aguas superficiales continentales. Ph. D. Thesis, Universidad Complutense de Madrid, Spain. Sihna, S.K., Holly Jr., Forrest, M., Dyer, M., Parris Lii, J.B., 1998. Threedimensional numerical modeling of thermal stratification in cooling ponds. International Water Research Engineering Conference Proceedings, vol. 2, Memphis, USA, 3–8 August, pp. 1044–1049. Wetang’ula, G.N., 2001. Ecological risk assessment of Nesjavellir Geothermal Power Plant wastewater disposal in lake Thingvallavatn in SW-Iceland. Geothermal Training Programme, Reports, 2001, The United Nations University, number 16 www.os.is/solofile/1908. Wu, J., 1969. Wind stress and surface roughness at air-sea interface. Journal of Geophysical Research 74 (2), 444–455. Wunsch, C., Ferrari, R., 2004. Vertical mixing, energy and the general circulation of the oceans. Annual Review of Fluid Mechanics 36, 218–314.