Experience of employment of computational models for water quality modelling

The purpose of the article is to develop the required and sufficient conditions under which numerical methods can be used for engineering calculations and for scientific research of hydrodynamic processes in solving practical problems related to predicting the spread of pollutants in water bodies and streams. The conducted studies consisted in comparing the results of laboratory experiments and mathematical modelling, in particular the distribution of heat in a stream with different temperature in water layers was studied. To check the adequacy of the proposed numerical models, calculations were performed and comparisons were made with the results of experimental data. The obtained results allowed to determine the boundaries of the qualitative difference in the flow behaviour for different numbers of Froude and Reynolds. The accuracy of the method was also studied. A number of additional requirements for numerical models were proposed in addition to approcsimation and stability, such as requirements of conservativeness (divergence), existence of trivial solutions on grids, possibility to calculate highly unsteady, quasi-stable, pulsating and stationary flows, requirement of invariance of linearized equations, as well as the requirement of a one-dimensional scheme to be a consequence of a twodimensional scheme. Distribution of velocities of wind currents using a three-dimensional and two-dimensional model was studied for a real object. A shallow-water bay of the Aral Sea was chosen as the object for the research. Comparison of the calculation results for both models showed that the flow velocity fields, as well as the distribution of pollutants in shallow waters, can be performed using a two-dimensional model. * Corresponding author: krutov@bk.ru © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). E3S Web of Conferences 97, 05030 (2019) https://doi.org/10.1051/e3sconf/20199705030 FORM-2019


Introduction
Increase in the volume of use of water resources and contemporaneous increase of the requirements for the quality of water withdrawn set a number of tasks the solution of which can impact on the efficiency of the facilities and environmental safety. The most important of the emerging tasks are:  determination of velocity fields in lakes and streams;  determination concentrations of pollutants, including temperature and salinity in lakes and streams;  prediction of bed flow deformations, roiling, and sedimentation;  determination the degree of eutrophication of water bodies;  development of engineering measures to regulate water quality both in water bodies and at water intakes;  forecasting of water quality in water bodies. The article considers options to cope with one of the mentioned above tasks as movement of concervative (nonpartitionable) pollutants in surface wates and consists of the following parts: basics on hydrodynamics equations used for mathematical modelling of hydraulic processes in water bodies, discussions of numerical models, major requirements for numerical models, comparison of the results of numerical calculations with experimental data, and conclusions.

Basic equation
In general, movement of pollutants in surface waters is an extremely complex combination of hydrochemical, hydrobiological and hydrodynamic processes. A complete picture of the quality of water cannot be described without taking into cjnsideration chemical and biochemical kinetics, sedimentation and diffusion of, for example, organic matter and a number of other physicochemical processes. In this paper, only hydraulic aspects of the movement of pollutants in surface water sources are considered and it is assumed that the interaction of chemical components occurs at a high rate, much higher than the propagation velocity of small perturbations (c=(gh)1/2+|U|, here h, |U| the depth and modulus of the flow velocity, respectively) and is described by first-order equations. Sufficient accuracy is required for the design and operational practice as well as high efficiency of the methods to apply at engineering practice. These conditions are the prerequsits to solve problems while performing mass calculations of investigated objects assessing their current condition and forecast possible developments. Currently, the most effective way of combining the mentioned conditions is mathematical modeling using the equations of hydrodynamics, which allows to identify the most important parameters necessary to obtain their numerical values for environmental protection practice. The solution of the general equations of hydrodynamics is a complex problem and there are still no effective algorithms for solving three-dimensional equations in full definitions. However, there are solutions to overvome these difficulties by introducing some hypotheses.
According to practice [1, 3 -19, 30 -32], for water badies whose horizontal dimensions are much greater than the depth, this can be done by introducing a large scale to consider the phenomenon of pollution movements. Depending on the degree of model upscaling, threedimensional baroclinic fluid equations, two-dimensional Saint-Venant equations, onedimensional equations and zero-dimensional equations (balance) could be considered. In this paper, to study the capabilities of numerical models that reflect physical processes in water bodies and streams the following system of equations was used [2]: here u i -projection of the current velocity vector onto axis x i , w -projection of the current velocity vector onto axis z,  -liquid density, g -component of the acceleration vector of gravity force,  -kinematic coefficient of viscosity, D -coefficient of vertical diffusion, S r -average substance concentration, q Sr -internal sources of a substance for the Newtonian fluid.
The essence of the method of finite differences consists in discretization of time and space variables (grid construction) and unknown functions (introduction of grid functions) from derivatives and integrals of these functions (construction of a difference scheme on the grid) [3 -19]. As a result, instead of a system of differential or integral equations at each point of the considered area, a system of algebraic equations in each discrete structure (cell) of the considered area will be obtained.
2 Theoretical preliminaries: Basic requirements for numerical models 1. Conservatism (divergent). The requirement of conservatism is widely reflected in literature [4 -7, 9 -20, 22 -30 ], including the one related to the establishment of boundary conditions [2]. 2. The scheme must have a trivial solution. On grids that do not change over time, the fulfillment of this requirement is not difficult. If, however, the grid on the real region changes with time, then this requirement leads to the necessity of introducing additional terms into the difference equations [3,14,36] 7. A one-dimensional scheme for straightened channels, taking into account the real change in the cross-flow area, would be a consequence of a two-dimensional one [1,9,11,19]. The following conditions required for vodelling of pollutions transfer [2]: 1. Divergences (the caracteristic to accurately retaining the grid analogy of the impuritymass).
2. Invariance of the difference equations with respect to the transformation S S C    , here S -pollution concentration, С -constant, 3. Symmetry. 4. It is necessary to comply with the independence of the stationary state from time.

Numerical models
Following [8] complete equations of motion (without the hypothesis of hydrostatics): ..
In the equations (4), it was assumed that the variation of the stresses along the horizontal coordinates is much smaller than along the vertical. Having made such an assumption about the flow velocities, the following was obtained from (4): . In order to consider the task in the vertical limits independent of x, and y a new coordinate system could be introduced: ; ; in which the flow area lies in the range from '0' to '1'. In these coordinates, the equations have almost the same form as in the old coordinates, with W replaced by W 0 : Henceforth , the sign "" over variables is implied.
In order to find fundamental approaches and cope with the task, let us assume that  and  are the functions which weakly depend on coordinates and time, and  is the known (set) function. In addition, let us neglect the convective vertical transport in comparison with the diffusion. These assumptions do not affect fundamentals of the proposed approach. Moreover, a large number of tasks can be solved with these assumptions [2,8]. Thus, under the assumptions made, we have the following system of equations: With boundary conditions: Obviously, conditions (6.1) -(6.3) are not sufficient. In fact, as reported in [21], this is the main problem when calculating a turbulent flow. This problem is referred in the literature as the 'enclosure problem'. One more condition is necessary. However, there are no physical prerequisites for its setting. This happened because by moving to the equations in 'stresses', we eliminated 'pressure', i.e. the additional condition is laid down in the original statement of the task. Indeed, integrating the first two equations (1), the missing condition could be obtained: here: With the introduction of equations averaged over the depth, the enclose problem is solved. In this case, it turns out to be more convenient instead of condition (6.3) to use the averaged equation of continuity: An algorithm for solving equations (5) for flows with high viscosity using Chebyshev polynomials was proposed in [2]. In [2], numerical experiments were carried out to investigate the accuracy of the method in the following way: for a given number of operations M=nN various ratios of the number of layers 'n' were taken and the number of polynomials 'N': 1. n=1, N=21 2. n=7, N=3 using the heat equation: with the initial condition . It was shown that the best option is to use a quadratic polynom in combination with the Galerkin method [2]. One of the serious problems encountered in the prediction of the spread of pollutants is the problem of adequacy of simulated results to actual processes. In this respect, the only criterion is the comparison of calculations with actual observations.

Comparison of the results of numerical calculations with experimental data
A large number of studies were devoted to the problem of adequacy of the models which describe hydrodynamic phenomena in water bodies, for example [10 -16, 26 -30,25-27].
To verify the adequacy of the proposed models, calculations were performed and comparisons were made with the results of model studies published in [34].
To analyse the effects of stratification on flow behaviour, velocity and density measurements were taken under different flow conditions. The techniques used for these measurements are explained below.
The essence of the experiment [34] was as follows: the rectangular channel consisted of a 50.0m long flume, whose first 0.5m were adapted to create a two layer stratified flow. A separate supply of two flows with different temperatures moving at different speeds was arranged by two inlets. The height of each inlet was 0.1 m (Fig.1).
After 0.5m the channel opened and the fluids from the two ducts started to mix. The depth of the flow was controlled by an outlet weir to obtain uniform flow with a flow depth of 0.2 m. These dimensions classify the channel as a narrow channel, with aspect ratio (BIH) of 1, and relative length (LIDh) of 100. T   Fig. 1. Scheme of the model.
The average temperature was measured continuously using a platinum resistance of the following dimensions: length -1.5 cm, width -0.3 cm, connected to secondary equipment. The average velocity was measured by micro-curent meter of 0.9 cm diameter every 2 cm from the bottom to the surface. The averaging time was 100 s. The mean square deviation of the temperature fluctuations was measured through each 1 cm of the profile using a heated wire sensor with dynamic feedback and integrator, which allowed to catch information at frequencies from 0.2 to 10 hertz. This seems sufficient for the experiment as these frequencies are the most important. The effect of the following independent parameters was studied: Experiments allowed to establish boundaries of the qualitative difference in the behavior of the flow at various Froude and Reynolds numbers. It was found that at Fr  1,6 there is a transition from a stable stratification to an unstable. Fig. 2 illustrates the comparison of the experimental data and the numerical implementation of the same model for two Reynolds (5*103 and 1*10 4 ) and Froude (0.9 and 5.0) numbers.
As a result of the experiment, vertical profiles of mean temperature (T), average horizontal speed (U) and mean-square deviation of the temperature fluctuation amplitude were identified. As can be seen from the figures, the velocity and temperature of water, determined by calculations, best of all coincide with the experimental results for case a), i.e. under Re = 5 * 10 3 and Fr = 0.9. In this case, the maximum deviation of the calculated water velocity on the surface from the one measured does not exceed 0.20 U/U 1 in the whole experimental area, and at the level y/h = 0.9 at a distance x/h = 30 it was 0.25 U/U 1 . In case b), i.e. under Re = 1 * 10 4 and Fr = 5.0, the maximum deviation of the calculated water velocity from the measured one on the water surface does not exceed 0.05 U/U 1 , and at the level of 0.8 y/h does not exceed 0.25 U/U 1. The maximum deviation of the calculated water velocity was fixed at a distance x/h = 30 at a depth of 0.9 y/h. At this point, the deviation was 0.34 U/U 1 . At a distance x / h = 100, the coincidence of the flow rates obtained in the experiment and as a result of the calculations turned out to be good in both a) and b). The proposed model better described the temperature distribution in the flow with Fr = 5.0. Most likely this is due to the fact that the model describes the turbulent flows fairly well. The problem of adequacy lies not only in comparison of the results of calculations with the data of experimental studies but to a large extent in the complexity of conducting both laboratory and full-scale experiments [29]. From this point of view, the results of the comparisons given in [9,2] show that models similar to (1) are in a good agreement between the results of measurements and calculations. A good agreement between the experimental data published in [27,28] and the numerical implementation suggests that the proposed model describes turbulent flows well and can be used for numerical modelling of objects with different-density flows. Comparison of acual wind borned current in the channel in physical model [34] and calculations were made to of determine the adequacy of the proposed model. The result of the comparison is presented in Figure 3.

Conclusions
The accumulated practical experience of applying finite-difference schemes for the solution of dynamic tasks made it possible to find the necessary conditions under which numerical methods can be used for engineering calculations and studies of hydrodynamic processes. To solve practical problems associated with predicting the muvement of pollutants in water bodies and streams, a grid-spectral method was developed that showed high efficiency in the case of stratified currents. This is especially justified for solving problems associated with the calculation of stratified flows. An algorithm based on the use of integro-interpolation methods was developed. The algorithm allows to carry out calculations with large time steps. Numerical studies using a three-dimensional model of flow and dispersal of pollutions in the investigated reservoir showed that, in general, the velocity fields at different depth levels are identical, i.e. the influence of three-dimensionality is insignificant. Therefore, it is advisable to carry out numerical modeling of the distribution of pollutions across the water area of reservoirs with similar characteristics using a two-dimensional (plan) model. Numerical simulation using a two-dimensional model gave a picture of the velocity and impurity fields over the area of the reservoir in the dynamics of development over a period of 5 years as well as the field of maximum values of pesticide -HCH.