Uncertainty in computations of the spread of warm water in a river – lessons from Environmental Impact Assessment case study

The present study aims at the evaluation of sources of uncertainty in modelling of heat transport in a river caused by the discharge coming from a cooling system of a designed gas-stem power plant. This study was a part of an Environmental Impact Assessment and was based on twodimensional modelling of temperature distribution in an actual river. The problems with the proper description of the computational domain, velocity field and hydraulic characteristics were considered in the work. An in-depth discussion on the methods of evaluation of the dispersion coefficients in the model comprising of all four components of the dispersion tensor was carried out. It was shown that in natural rivers all components of a dispersion tensor should be taken into account to qualitatively reflect the proper shape of temperature distributions. The results considerably depend on the 2-D velocity field as well as hydraulic and morphometric characteristics of the flow. Numerical methods and their influence on the final results of computations were also discussed. All computations were based upon a real case study performed in Vistula River in Poland.


Introduction
A prerequisite for the construction of a new industrial plant is an Environmental Impact Assessment (EIA), understood as a formal process used to predict the environmental consequences of any development project.In case of gas-stem power plants, as part of the EIA, one needs to evaluate the impacts caused by heated water discharged into a river.It is rather obvious that the credibility of the computations performed within EIA depends strongly on the available data.In practice, however, such data is often of limited reliability, yet quantitative results from relevant models are expected.In other words, time constraints and lack of information require that the EIA must rely exclusively on expert skills and opinions.Moreover, planning of additional measuring campaigns during EIA may become impossible within the given time.Therefore, the discussion of the usefulness of the obtained results is of crucial importance.The authors aim to share their experience in the use of an up-to-date model on the spread of a warm water jet in a river in the light of scarcity of proper data.It is quite a common practice that professionals gain personal experience, particularly with respect to the models developed by themselves.Usually they have good knowledge on how their parameters and boundary conditions should be assigned to yield satisfactory results, and at the same time they realise that their model may fail under some circumstances.Evaluation of such circumstances constitute, in fact, an empirical estimate of the uncertainty of the simulated results.Beven (2007a) claims that in non-ideal cases (i.e., nearly all real applications) non-statistical (epistemic) uncertainties may dominate.They are e.g., bias and nonstationarity in input errors, model structural errors and commensurability errors (where a variable or parameter in a model is different to an equivalent quantity that can be measured in the field).In practice, it is almost impossible to separate different sources of aleatory (statistical) and epistemic uncertainties unless very strong assumptions are made (Beven, 2007b;Beven et al., 2011).Nevertheless, assessing the sources of uncertainty with respect to the processes under study is crucial for the analysis (e.g.Catari et al., 2011;Hughes et al., 2011).Only a few sources of uncertainty, but those which significantly influence the results, will be studied herein.It is by no means a formal uncertainty analysis -even simple Published by Copernicus Publications on behalf of the European Geosciences Union.
goodness-of-fit criteria cannot be applied because no observed values are at our disposal to compare with the predicted dependent variables.In uncertainty analysis there are obviously other possibilities, for example, with use of Monte Carlo-type simulations that are capable of predicting model output variability from a large number of deterministic model runs, with input variables sampled from appropriate distributions (e.g.Kochanek and Tynan, 2010;Scharffenberg and Kavvas, 2011;Shen et al., 2012).The question arises, however, as to which distributions would be appropriate for the particular case considered.Moreover, the computational expense of multiple model runs would definitely be too high.
One source of uncertainty is insufficient knowledge of the analysed phenomena.In addition, we have to consider errors introduced by the models used in the calculations, which always simplify the described phenomena.For various reasons these models cannot take into account all the variables controlling the phenomena, and usually simplifications are made to make the problem practically solvable.Possible numerical errors caused by the applied numerical algorithms have to be considered as well.
The analysis of extreme cases, in which the natural environment is endangered most, is very often considered in EIA as the worst-case scenario, in other words, as the worst performance among all scenarios that could be generated.
The present study is based on a case study aimed at building scenarios of the spread of heated water discharged from a designed gas-stem power plant on the Vistula River below Włocławek town (Kalinowska et al., 2012).Such heated water constitutes an environmental problem in many situations and, therefore, is often called thermal pollution.It is usually a side effect of power plants operations, chemical industries and other hydraulic engineering facilities using water from rivers and open channels in the cooling process.We encountered a number of problems significantly influencing the credibility of the results, which motivated us to share this experience with the readers.In the study, a two-dimensional temperature field resulting from a warm water discharge was sought.

Mathematical model of heat transport
When developing models there is always a trade-off between retaining greater detail of the processes under study and deriving tractable equations.Inevitably, some details are simplified in, or even excluded from, the models used in practice.Although, in general, the temperature distribution in rivers should be described by the three-dimensional (3-D) heat transport equation, its reductions to two (2-D) or even one dimension (1-D) are considered in practice.The main obstacles to use a 3-D approach is the lack of knowledge of the realistic detailed 3-D velocity field and high computational expense.In the so-called mid-field region (stretching down the river until complete lateral mixing occurs) the depth-averaged heat transfer models are relevant (Rodi et al., 1981;Seo et al., 2010;Szymkiewicz, 2010): where: t -time, x = (x, y) -position vector, T (x, t) -depthaveraged water temperature, h(x) -local river depth, v(x)depth-averaged velocity vector, D(x) -heat dispersion tensor, q(x, t) -source function describing additional heating or cooling processes.This equation represents the temperature change in time due to heat flows: that is heat carried with average velocity, carried by velocity fluctuations (turbulent heat conduction), transmitted by the deviation of velocity and temperature across the depth (dispersion of heat) and heat flows to and from sinks and sources, respectively.The dispersion tensor that appears in the equation due to the depth-averaging represents an additional, significant transport mechanism, which is not a physical process, but only a consequence of averaging of the equation.By analogy with the turbulent heat exchange, it is assumed that dispersion flow is proportional to the average temperature gradient (Rowiński, 2002;Rutherford, 1994).This proportionality is determined by the so-called dispersion coefficients (D xx , D xy , D yx , D yy ).It is worth noting that similar molecular and turbulent heat conduction coefficients are several orders of magnitude smaller than the dispersion coefficients.Usually the molecular coefficients are omitted and the turbulent ones are either omitted or included in the dispersion tensor.
Depth-averaging results in model simplifications.Water in the rivers is usually well-mixed and the temperature is nearly uniform from surface to bottom (Allan, 1995), but nevertheless we should always remember that some variations caused by external factors (e.g., surface water temperature exchange with the atmosphere, groundwater, etc.) may occur.Sources and sinks of heat energy may play an important role in Eq. ( 1), but in the case considered herein we lacked information on these components.They embrace the exchange with the atmosphere and with the groundwater through e.g., shortwave and longwave radiation, latent heat transport, ground heat flux, etc.To some extent we may assume that these additional terms affect both the heated water and the ambient water in more or less the same way, and since we are interested in the difference between the temperatures in the river and in the warm water jet, the source term may be neglected in the computations.But this simplification is definitely a potential source of error that may arise in computations.On the other hand, inclusion of this term with guessed values could introduce more serious errors.

Solution to 2-D heat transport equation
To solve the heat transport equation described in the previous section, we need to know: the geometry of the river, the two-dimensional velocity field, the boundary and initial conditions, the full dispersion tensor, and if we want to include any additional sources (e.g., temperature exchange with the atmosphere), substantial additional information is required (like meteorological data).We may encounter serious problems in the process of data collection and acquisition which frequently lead to simplifications affecting, to some extent, the final solutions.When solving typical academic problems (especially in laboratory), one usually tries to gather all necessary data and consider all possible processes which affect the solution.In real situations when dealing with EIA, such an approach is often practically impossible.Measurements of potentially necessary data are usually limited due to time constraints, huge costs and many technical restrictions.Recently Piotrowski et al. (2011) have analysed methods for analysis of the fate of pollutants over long distances (far field) for applications in ungauged rivers, i.e., where little hydraulic or morphometric data are available.In the present study, a different approach is proposed and it is limited to the mid-field only.However, the obtained solutions have errors that we should be aware of.We need to have in mind what aspect is of crucial importance for the decision-maker requesting the EIA and, in this regard, it may be sufficient to show the possible range of the values of interest (and not the precise numbers).Knowing the possible minimum, maximum and mean values we may be able to prepare realistic scenarios (including the most significant extreme cases).
Below are discussed some problems that may be encountered in the preparation of a scenario of the spread of heated water discharge from a power plant.Various aspects from an actual EIA study (selected for the purpose of illustration) are considered herein.

Problem considered
This study is based upon the computations of the spread of the heated water discharged from a designed gas-stem power plant located near Włocławek town on the Vistula River.The variability of water temperature depending on four locations of differently designed exit pipes was of particular interest within the study.The water is supposed to be discharged with a constant flow rate of 14 m 3 s −1 and a temperature that is 7 • C higher than the temperature of the ambient water.The adopted water flow in the river was the averaged low-flow: Q = 334 m 3 s −1 .This is a condition with almost minimum dilution of the heated water discharge and, hence, the largest temperatures.The detailed description of the study area, the discussion of the computations and final results may be found in Kalinowska et al. (2012).

Geometry
The first problem is encountered already at the level of the definition of the computational domain.This is a well-known problem in hydraulic computations, but the problem is not often revealed by modellers when they publicise their results.It is worth noting that characterisation of the computational domain for natural rivers is a type of art even in 1-D situations when the 1-D velocity field is the main concern of the modeller (see e.g.Rowiński et al., 2005b).
A particular problem in this study was related to an insufficient number of measured cross-sections, which required special interpolation procedures to be used to generate consecutive transverse depth profiles for the model.Assuming that cross-sections were carefully chosen accounting for all the critical points (like river bends, dunes, river narrowness, etc.), which obviously is not always the case, the bathymetry and the computational grid can be described mathematically.In the considered case 21 cross-sections of the Vistula River (between the 690.250 and 718.200 km of the Vistula River) were available.The problem with them was that they were measured a relatively long time ago, namely in 1994.The data obviously does not truly represent today's situation.An example of the variation in the bed elevation during an earlier period (between 1971 and 1994) for a selected cross-section is shown in Fig. 1.This change is significant, being as large as 2 m, and one may expect that equally significant changes could have occurred after 1994 as well.
Since some preliminary computations showed that the analysed stretch of the river may be reduced (transverse mixing is faster than initially expected), only 9 cross-sections were used for the final calculations.The average distance between the cross-sections turned out to be as large as approximately one kilometre and obviously some of the variability in bed elevations was omitted.
It is quite obvious that the number of the measured crosssections was too small for describing the shape of the river channel adequately and, therefore, it was necessary to generate a topographic map and to create (by interpolating procedures) additional cross-sections before the final bed profile and computational grid were produced.The two-dimensional computational grid and the bed profile were prepared using the CCHE MESH generator developed by NCCHE -National Center for Computational Hydroscience and Engineering (Zhang, 2005a).This is a user-friendly mesh generator for generating structured quadrilateral meshes based on the bed topography and the bed elevation data.The grid was then used to calculate the velocity field by means of the CCHE2D model (Zhang, 2005b), also developed by NCCHE.A different type of grid was required by the River Mixing model (RivMix), developed by the authors, used to compute the temperature distribution (Kalinowska and Rowiński, 2008).This is often the case when two different computational models are used.Then all necessary data have to be transformed from the computational domain of one model to the domain of the other one.This procedure should be performed carefully, because this is another stage where additional errors may be introduced.Proper interpolation procedures -Delaunay interpolation (de Berg et al., 2008) -were used to transform the data in the analysed case from one computational grid to the other.In this way, a number of grids with different grid steps have been prepared.For the final computations and analysis presented in the paper a grid with step sizes of x = y = 10 m was chosen.Such step sizes, together with an appropriate choice of numerical method, gave us a fast and sufficiently accurate solution in the considered situation.Figure 2 presents the resulting interpolated water depth on the chosen grid.
The water depths were computed for the selected discharge and water level elevation.It should be noted that any change in these parameters would definitely affect the solutions.In this case study, it was impossible to perform calculations for different ranges of possible water levels and discharges and, therefore, it was assumed to be sufficient to carry out the calculations for a reasonably low water level to capture the environmentally most severe situation.
It is also important to note the limitations of the twodimensional model.Since 2-D depth-averaged models can be only used after the complete vertical mixing occurs, we cannot interpret any results in the initial stage in the direct neighbourhood of the discharge (so-called near-field zone).According to the procedure proposed in Jirka and Weitbrecht (2005), this initial distance in the considered case reaches about 75 m from the discharge site.

Velocity field
The knowledge of a relatively accurate two-dimensional velocity field is a prerequisite to solve the two-dimensional heat transport equation (Eq.1).In reality we usually do not have enough measurements, especially when dealing with a 2-D case; therefore, it is necessary to model the velocity profile.There are numerous software packages available, but nevertheless it is essential to bear in mind that the solution of the momentum equations is not straightforward.It is necessary to know what kind of equations are solved in the applied software and what kind of numerical algorithms are being used.Then we may have some expectations on the resulting errors.In the considered case the CCHE2D model -twodimensional depth-averaged, unsteady turbulent open channel flow was used (Altinakar et al., 2005;Jia and Wang, 2001;Ye and McCorquodale, 1997;Zhang, 2005b).The model is based on the depth-averaged Navier-Stokes equations.The details of the velocity field are given in Kalinowska et al. (2012) and will not be discussed herein.However, one needs to remember that an accurate evaluation of the velocity field is crucial for further computations and it is burdened with all the problems that turbulence modelling bring, such as choice of the so-called closure hypotheses, etc.The resulting velocity magnitude distribution is presented in Fig. 3.

Initial and boundary conditions
The nature of differential equations means that one needs to know, a priori the initial and boundary conditions for the problem to be solved.In this study, an average natural river temperature of the considered area was taken as the initial condition.It was assumed that the depth-averaged 2-D river temperature was the same within the entire domain area.We do realise that in reality the river temperature may change in space or in time due to day-night or seasonal changes.Taking into account those changes would require yet more experimental data.Since we are interested in the water temperature distribution after the release of the heated water (more precisely: increase of the water temperature values), calculations may be performed with regards to the relative temperature ( T ) and may readily be scaled to the river temperature at a given time of a day or a season.
In the considered case the continuous discharge of 14 m 3 s −1 of water heated by 7 • C in relation to ambient water was assumed in the computations.It was assumed that the river bank had a temperature equal to the initial temperature of ambient water, and the boundary conditions were defined in a way such that the river water was not allowed to heat the bank, but the bank could possibly cool down the water.Such assumption was caused by the lack of information on the bank temperature, but could locally influence the results.It may particularly be important when the source of the heated water is located on the river bank and the difference between the bank and the source temperature is large.
There are numerous ways of discharging the heated water to the river and finding those that minimise the environmental impact is a key element of the EIA.Several options were analysed in the study, varying according to the location and method of discharge: a continuous point-like discharge through a single nozzle or a continuous distributed discharge through a series of uniformly spaced nozzles along a straight exit pipe.Four of those variants -two point-like at different locations and two distributed with exit pipe lengths of 14 and 28 m -were presented in Kalinowska et al. (2012).
The numerical implementation of the discharge is again only a poor approximation of the real conditions.Solving numerically the heat transport equation, the continuous solution domain has to be replaced with a discrete domain and the line segment pipe within the considered grid is represented by discrete points in the grid nodes located along the segment.In case of the grid with steps of 10 m, the 14 m line segment is represented by just 2 points (located near the beginning and the end of the segment) and the 28 m line segment is represented by 3 points (near the beginning, the centre and the end of the segment).Such approximations are obviously a major simplification.An additional serious simplification is made in the process of calculation of the effective temperature (T E ) at a relevant discharge site.Since the single grid cell of dimension of x by y is the smallest fragment of the river that we can consider, it was assumed that in such single cell surrounding the point of discharge the heated water mixed immediately with the river water.Taking into account the volume of the river water and the heated water in the considered cell we calculated the effective temperature after the discharge based on a simple expression: where: T W -the temperature of the ambient water, T Z -the temperature of the discharged heated water, Q W -river flow at the source, Q Z -heated water flow at the source, Qthe resultant (total, i.e., Q W + Q Z ) flow discharge.The river flow discharge at a given point at the source was determined by the relationship: where: h Z -the water depth at the source, v Z -water velocity at the source, x -grid spacing.Simple calculations show that the effective temperature obviously depends strongly on the assumed grid spacing, and the effective temperature naturally tends to 7 • C when the grid size approaches zero (see Table 1).
The above assumptions involve no loss of generality -they influence the results of the computations mainly in the close surroundings of the discharge location, an area that cannot be properly interpreted by the 2-D model anyway because the processes here require a 3-D treatment.

Dispersion coefficients
Dispersion coefficients controlling the rate of mixing are essential to solve Eq. (1) and at the same time their determination is extremely complex.In the general case using Cartesian coordinates the dispersion coefficients form a nondiagonal dispersion tensor with four dispersion coefficients: Those four dispersion coefficients should be computed on the basis of the so-called longitudinal D L and transverse D T dispersion coefficients.It is crucial to compute the dispersion tensor D in the proper way, otherwise unrealistic results may easily be obtained.The authors in their earlier studies (Kalinowska and Rowiński, 2008;Rowiński and Kalinowska, 2006) presented how erroneous ways of simplifying the treatment of the dispersion tensor often met in literature, result in misleading concentration distributions revealed both in their shapes as well as the concentration values.In brief, the proper way to obtain the full tensor D is rotation of a diagonal tensor D D containing the D L and D T coefficients: where: R = cos α − sin α sin α cos α -rotation matrix, α -angle between the flow direction and x-axis, D D = D L 0 0 D T .There are several erroneous ways of simplifying the tensor, but they usually boil down to somehow neglecting the off-diagonal components: D xy , D yx .Figure 4 presents the distribution of the relative temperature T (the difference between the temperature of ambient water and the actual river temperature) in case of the point-like continuous discharge in the middle of the channel along the cross-sections located ca. 250 and 500 m from the discharge point.The results were obtained for four different methods of computation of the dispersion tensor D: using the proper method defined by Eq. ( 5) and marked as I; using simplified methods marked as II (the off-diagonal elements of the dispersion tensor are omitted), III (dispersion coefficients D L and D T are treated as a vector), and IV (the diagonal elements of dispersion tensor D xx and D yy are simply replaced by D L and D T , and the offdiagonal elements are 0).The difference between the results can be easily observed.The increase of the temperature in the considered case in the middle of the channel is much bigger (almost 3 times) when using the proper way of calculation of the dispersion tensor (I), which could be very important in case of an EIA.Also the mid-field zone is much larger in this case.In case of the simplified method IV, the cloud of thermal pollution is pushed to the left river bank.Detailed definition and analysis of the simplified methods can be found in Rowiński and Kalinowska (2006).
Longitudinal and transverse dispersion coefficients depend on many factors related to the geometry of the channel, dynamics and turbulence of the flow.Due to their significance and difficulties associated with their determination these coefficients are subject of many discussions in literature, but there are still a lot of questions and misconceptions concerning their values.The best way of estimating the dispersion coefficients for the actual river is a tracer experiment, but usually in practical applications -since such an experiment is expensive and time consuming -impossible to be carried out.Moreover, the practice shows that tracer tests are feasible in one-dimensional situations only and they allow for the determination of just the longitudinal dispersion coefficients (Deng et al., 2001(Deng et al., , 2002;;Guymer, 1998;Kumar and Dalal, 2010;Rowiński et al., 2008;Sukhodolov et al., 1998).
In cases when such tracer tests are not available, reliable estimation of dispersion coefficients becomes extremely difficult and can be a source of large uncertainty.Universal formulae for these coefficients have been sought in many research centres and they are usually related to the known hydraulic parameters, such as averaged depth H , width B and velocity U , shear velocity U * and the channel sinuosity S or the radius of curvature R of the considered channel.The derived formulae are of rather limited universality and often the formulae working for one channel do not hold for others.Another problem for the modellers is that many papers proposing or applying such formulae do not discuss their limitations or ways of their determination.According to selected review articles presenting expressions for transverse dispersion coefficients (e.g.Deng et al., 2001;Jeon et al., 2007;Seo and Baek, 2008), the values of D T for the considered area of Vistula River are plotted in Fig. 5. Bulk hydraulic parameters were used in all formulae and one can readily note differences that these various formulae produce.These differences are caused by a number of factors.In principle various formulae are constructed for different types of rivers and flow conditions and they take into account different processes and different variables.For example, the formula proposed by Fischer (1967) considers the transverse turbulent diffusion only.The formula was built for an artificial uniform and straight channel with constant depth.Other formulae take into account both the transverse turbulent diffusion and the transverse dispersion.In the Rutherford formula (Rutherford, 1994) it is assumed that the dispersion effect is significantly larger than turbulent diffusion which is omitted.Note also other problems that should be tackled when using these formulae.Some of them include, for example, geometrical parameters such as sinuosity index or curvature radius of the channel.Those values are easily measurable in laboratory conditions, but in natural rivers their determination may become extremely difficult and not unique.Meanders and bends may change significantly from one section to another and then the values of the dispersion coefficients should also change.Single values assumed for the whole river reach may introduce serious errors.Moreover, in practice the division of the river reach into sections with constant sinuosity may become intractable.The situation turns out to be not much simpler when the mean river width or depth are taken into account.In the considered case herein the average width of the channel is about 400 m.The values of D T coefficients obtained with use of the average river width are presented in Fig. 5 in green.The geometry of the considered Vistula River reach is such that at some cross-sections the river narrows to 300 m or even to 200 m.Then the values of D T change and these changes are revealed through the width dependence of the various formulae (see Fig. 5 blue and red bars).
In case of the longitudinal dispersion coefficient D L , the number of formulae presented in literature is extremely large.Wallis and Manson (2004) review and discuss many of them.Numerous and successful attempts pertain also to other methods like, for example, artificial intelligence methods allowing to estimate this coefficient on the basis of the known hydraulic parameters (see e.g.Kashefipour et al., 2002;Piotrowski et al., 2011Piotrowski et al., , 2006;;Rowiński et al., 2005c;Tayfur and Singh, 2005).
It turns out that the differences between the values of the obtained D L , with use of different methods or expressions, are larger than in case of the coefficient D T .The differences in the obtained D L values by means of various methods are often one or more orders of magnitude.There is also another crucial problem related to the basic definition of the D L coefficient.The longitudinal dispersion coefficient appearing in 2-D equations is not the same entity as the longitudinal dispersion coefficient in 1-D equations and this fact is very often forgotten in scientific considerations.The values of dispersion coefficients relevant in 1-D situations are very often adopted to 2-D models which on one hand is often done unconsciously, and on the other hand is a practical approach since usually only values for the 1-D approach are available.This introduces additional and eventually quite serious errors.Let us illustrate the above problem.In the case study considered in this paper there was no experimental data that could be used to calibrate the model to obtain the proper values of the D L and D T coefficients.As expected, different relationships for dispersion coefficient resulted in significantly different values of D L and D T .In the study the most general relationship for dispersion coefficients (Czernuszenko, 1990;Sawicki, 2003) was used: where: D -longitudinal or transverse dispersion coefficient, u * -bed shear velocity, a -dimensionless parameter that theoretically may assume values from a relatively large range.
Taking into account the reasonable maximum and minimum values of the parameter a (based on the experience gained in similar rivers) we were able to analyse the environmentally most severe cases, which is actually the main task of an EIA.The most probable value of a was also analysed and presented in the study.The choice of a is not free from subjectiveness.
The value of the parameter a is not the same for longitudinal and transverse dispersion coefficients.Usually the ranges of a for D L and D T are the following (e.g.Rutherford, 1994;Sawicki, 2003): 30 ≤ D L h u * ≤ 3000; 0.15 ≤ D T h u * ≤ 0.9. (7) Note the relatively large range of the parameter a.In principle the range for the longitudinal dispersion coefficient (Eq.7) is provided for the 1-D case, but this range and especially the lower range, coincide with the values of a for the 2-D case.Since no information is available on the range in the 2-D case, the whole 1-D range is presented to cover a variety of admissible situations.When detailed geometric and bathymetric data is provided (which is rarely the case), it might be possible to obtain ranges of dispersion coefficients for 2-D models (Piasecki and Katopodes, 1999).Unfortunately, in the considered case (which is a frequent situation) the provided data is extremely scarce and we may work only with anticipated values.On top of that, any averaging of the value of h -either within the entire river reach or even between two given cross-sections may introduce a substantial error.Another problem of quite fundamental nature is the determination of the bed shear velocity (see Rowiński et al., 2005a).
In this study the bed shear velocity has been calculated based on the following basic expression: where: τ -total shear stress, ρ -water density, the calculation assumes ρ = 1000 kg m −3 .The distribution of the total bed shear stress was calculated with use of the CCHE2D model based on the turbulent Reynolds stresses.Since we do not know the value of the parameter a for the considered case in the present study, we conducted a series of simulations for different values of a.The variant with the point-like continuous warm water discharge of 14 m 3 s −1 , located in the middle of the channel has been chosen to illustrate the problem.Figure 6 presents the 2-D temperature fields for three selected values of parameter a used to calculate the longitudinal dispersion coefficient D L (a = 100, a = 500, a = 1000).Plots in the left column of Fig. 6 present the temperature for the whole area under consideration whereas plots in the right column zoom in on the area surrounding the point of release.Note that the plots present the relative temperature T -the difference between the actual river temperature and the temperature of ambient water.Figure 7 allows for accurate quantitative analysis of the differences between the distributions of water temperature calculated for various parameters a at the cross-sections located, respectively, at 100, 250 and 500 m from the point of discharge.For the convenience of the reader, these crosssections are also marked in Fig. 6 (right column).Note that in case of large values of the parameter a the temperature cloud slightly propagates also upstream.
Visible differences in the results reveal the importance of the selection of the dispersion coefficient.Therefore, potentially extreme values of a are considered to discuss the least favourable variant that may occur in the river reach.In comparison to longitudinal dispersion the selection of the transverse dispersion coefficient in the analysed case does not influence the results significantly and, therefore, computations of different variants are not presented.
Further, similarly to previous work (Kalinowska et al., 2012) the results of computations are presented for the   longitudinal and transverse dispersion coefficients calculated for a equal to 500 and 0.6, respectively, using the averaged values of the water depth (1.53 m) and the bed shear velocity (0.045 m s −1 ).In this case, the dispersion coefficients values are: D L = 34.425m 2 s −1 and D T = 0.041 m 2 s −1 .In case of using local values of the water depth and the bed shear velocity the distributions of dispersion coefficients are not uniform and are presented in Fig. 8.Note the differences in the solutions obtained with use of the averaged and local values of the depth and the bed shear velocities (see Fig. 9).The difference in results in the selected case is not significant.The maximum difference between the results (bottom chart)

Numerical solution
An important issue in the discussion of the model results is the mathematical consistency between continuum and discrete variants of partial differential equations used to represent the heat transport in a river.In academic studies an important step would be to compare numerical solutions obtained from calculations performed on successively refined meshes or grids to a reference solution.In the case considered herein, the exact solution of the advection-diffusion equations is not known to define such reference.Therefore, geometrically simpler problems are considered to assess the numerical solutions.In the evaluation of the numerical methods applied the key problems to be assessed are (Kalinowska andRowiński, 2004, 2008): numerical diffusion and numerical dispersion, instability and computational requirements.The first one -numerical diffusion -results in faster spreading of thermal pollution.It is a well-known effect, but since it may be difficult to observe in real applications, usage of numerical schemes that may generate such errors has to be made with care.Numerical dispersion errors may cause non-physical oscillations, but the effect is easier to notice in the obtained results.The details of the errors arising from the use of various numerical methods applied to Eq. ( 1) are given in Kalinowska and Rowiński (2007).In that study a two-dimensional model of the spread of passive pollutants in surface water -RivMix (which offers a choice of four different numerical schemes) was used.RivMix was also used in the current study, to model the thermal pollution, and an accurate and fast Alternative Direction Implicit (ADI) method was selected to solve Eq. (1).Following some preliminary numerical tests, appropriate time step ( t = 1 s) and spatial step ( x = y = 10 m) sizes were chosen to ensure a sufficiently fast and detailed solution.But before the time and spatial steps were chosen several numerical tests had been performed.The influence of the increasing model resolution on the final results is a subject of many studies in environmental fluid mechanics problems (see e.g.Ziemiański et al., 2011).

Conclusions
Mathematical models are suitable and powerful tools in Environmental Impact Assessment whenever influence of new hydraulic constructions on the natural environment have to be considered.Those tools cannot of course replace relevant observations and measurements and in the light of lacking data and impossibility of model calibration, understanding of various sources of uncertainties are crucial for the assessment of model credibility.This study was based upon an investigation of the spread of heated water discharged into a river from a designed gas-stem power plant.It was analysed how model simplifications may influence the final results.Particular attention was paid to methods of evaluation of dispersion coefficients.It was shown that in natural rivers all components of a dispersion tensor should be taken into account to qualitatively reflect the proper shape of temperature distributions.Omission of the off-diagonal dispersion tensor components may artificially lead to results which are better from the EIA point of view, i.e., the heated water would cool down more rapidly than in reality.Since the off-diagonal tensor components are obtained from the longitudinal and transverse dispersion coefficients, their determination is crucial for obtaining accurate results.Approximate values of dispersion coefficients for Vistula River could be calculated based on one of numerous phenomenological formulae, but only insitu tracer tests (which are logistically extremely difficult and expensive) could eventually provide credible results.Taking into account possible extreme values of these coefficients we are able to assess the most severe scenarios of the spread of warm water.The situation in which the worst scenario would not meet acceptable environmental standards might lead to carrying out measurements and/or tracer tests to determine the dispersion coefficients for the given river reach.The results depend considerably on the 2-D velocity field, hydraulic and morphometric characteristics of the flow, particularly the bed shear stresses.One should also be aware of the choice of numerical method that might introduce nonphysical phenomena in the final results.

Fig. 1 .
Fig. 1.Selected cross-section of Vistula River with marked bed profiles in years 1971 and 1994.Significant changes at the bed profiles are visible.

Fig. 2 .
Fig. 2. The water depth for the considered reach of the Vistula River, computed for the averaged low-flow discharge.

Fig. 3 .
Fig. 3.The velocity magnitude for the considered reach of the Vistula River, computed for the averaged low-flow discharge.

Table 1 .
The effective temperature in the single cell at the source in case of the point-like continuous discharge located at point Z = (1850, 775) for different sizes of grid cell.x= y [m] T E [ • C

Fig. 4 .
Fig. 4. Distribution of temperature ( T ) in case of the point-like continuous discharge in the middle of the channel along the crosssections located about 250 and 500 m from the discharge point: I -with the proper method of the dispersion tensor computation; II -with simplified method in which the off-diagonal elements of dispersion tensor D xy and D yx are omitted; III -with simplified method in which dispersion coefficients D L and D T are treated as a vector; IV -with simplified method in which the diagonal elements of dispersion tensor D xx and D yy are simply replaced by D L and D T , the off-diagonal elements are treated as 0.

Fig. 5 .
Fig. 5. Transverse dispersion coefficient for the considered reach of the Vistula River calculated using several formulae (taking into account different hydraulic parameters).For the convenience of the reader, the lower part of the chart is also presented in a different scale.b -dimensionless parameter; lab and field -denotes a test site (experiment performed in a laboratory channel or in a river).

Fig. 7 .
Fig. 7. Temperature distribution in case of continuous discharge of 14 m 3 s −1 of warm water at point Z 1 = (1850, 800) for different values of dimensionless coefficient across the cross-sections located at 100, 250 and 500 m from the discharge.
M. B. Kalinowska and P. M. Rowiński: Spread of warm water in a river

Fig. 8 .
Fig. 8.The longitudinal DL and transverse DT dispersion coefficients for the considered reach of the Vistula River, computed for the local values of river depth and bed shear velocity.

Fig. 8 .
Fig. 8.The longitudinal D L and transverse D T dispersion coefficients for the considered reach of the Vistula River, computed for the local values of river depth and bed shear velocity.

Fig. 9 .
Fig. 9. Two-dimensional temperature distribution results for continuous discharge of 14 m 3 s −1 of heated water along a straight exit pipe, 14 m long, located near the left bank of the river in case of using averaged (left-top chart) and local values of depth and shear velocity (right-top chart) for computing dispersion coefficients and difference between them (bottom chart).The beginning and end of the exit pipe are located, respectively, at points: Z 3B = (1730, 690), Z 3E = (1740, 700).