The program packages of realization of mathematical models of pollutants transport in rivers and for identification of river water excessive pollution sources located between two controlled cross-sections of the river will be considered and demonstrated. The software has been developed by the authors on the basis of mathematical models of pollutant transport in the rivers and statistical hypotheses test­ing methods. The identification al­go­rithms were elaborated with the supposition that the pollution sources discharge dif­fe­rent compositions of pollutants or (at the identical com­po­sition) different propor­tions of pollutants into the rivers. One-, two-, and three-dimensional advection-diffusion mathematical models of river water quality formation both under classical and new, original boundary conditions are realized in the package. New finite-difference schemes of calculation have been developed and the known ones have been improved for these mathematical models. At the same time, a number of important problems which provide practical realization, high accuracy and short time of obtaining the solution by computer have been solved. Classical and new constrained Bayesian methods of hypotheses testing for identification of river water excessive pollution sources are realized in the appropriate software. The packages are designed as a up-to-date convenient, reliable tools for specialists of various areas of knowledge such as ecology, hydrology, building, agriculture, biology, ichthyology and so on. They allow us to calculate pollutant concentrations at any point of the river depending on the quan­­tity and the conditions of discharging from several pollution sources and to identify river water excessive pollution sources when such necessity arise.

Keywords: software, pollutants transport, identification, excessive pollution sources, mathematical models

For solving the problems of analysis and control of the quality of the environmental objects, it is necessary to process big volume of measurement information about physical, chemical and biological parameters of these objects. To process big data with guaranteed quality in acceptable period of time is possible by means of wide application of mathematical methods and computers. For collecting, storage and processing data of the environment, there are developed the automated systems of the quality control of the environmental objects and universal software packages.1 Among them the most widespread are environmental water and air pollution levels control systems. Among the most topical problems of control of environmental water quality, the task of modeling of polluting substances propagation in water objects should be emphasized.

Transport and diffusion equation

The diffusion of polluting substances in rivers is mostly described by the three-dimensional equation of turbulent diffusion of non-conservative substances:^1-4

$\frac{\partial }{\partial t}\Phi (t,r)-\left(\nabla \cdot K(r)\cdot \nabla \right)\Phi (t,r)+\left(\nu (r)\nabla \right)\Phi (t,r)+\zeta (r)\Phi (t,r)=f(t,r)$ ,

Where $\Phi (t,r)$  is the time-averaged concentration of non-conservative pollutant; t is the time;  $r=[x,y,z]$  is the radius vector; it’s components  x,y,z are the spatial coordinates (the axis x  is horizontal and its direction coincides with the direction of averaged current of the flow, the axis  $y$  is perpendicular to the free surface and is directed downwards; the axis  $z$  is directed across the flow);  $K(r)$ is the tensor of turbulent diffusion,  $v(r)$  is the vector of time-averaged speed of the river flow;  $\zeta (r)$  is the coefficient of non-conservatism of pollutants,  $f(t,r)$  is the power of polluting sources.

For the engineering practice very often the two-dimensional or one-dimensional models are considered instead of the three-dimensional one.^2,5 If the non-uniformity of distribution of concentrations of pollutants on the depth of the watercourse is not taken into account, then it is possible to obtain the two-dimensional turbulent diffusion equation. The equation of one-dimensional turbulent diffusion is applied when the distribution of the concentration of the pollutant across the stream is homogeneous. It is also used when average indicators are used to represent pollution across the river.¹

Initial and boundary conditions

The unknown function  $\Phi (t,r)$  is defined at  $t>0$  and  $r\in G$ , where  $G$  is some region with the boundary  $\partial G$ , or  $r=x\in [0,L]$  (in one-dimensional model).

The additional conditions are specified in the form of

$\Phi (0,r)=S_0,\Phi (t,r){|}_{x=0}=\sigma$

( $S_0$ ,  $\sigma =\text{const}$ ). The boundary conditions in the lower extremity of the river section may be classical or non-classical. The classical conditions (condition of complete mixing) look like

$\frac{\partial }{\partial x}\Phi (t,r){|}_{x=L}=0$ ; (condition of full mixing)

non-classical conditions –

$\Phi (t,r){|}_{x=L}=q·\Phi (t,r){|}_{x=L-l}$ . (not local boundary condition)

Here  $q$  is the coefficient of self-cleaning of the river on the considered section,  $l$  is the length of the section.^5,6

At using of two- or three-dimensional models, the following boundary conditions on the other part of the line or on the surface  $\partial G$ (Neumann condition), are set also:

$(\nu ·\nabla )\Phi (t,r){|}_{r\in \partial G}=0$ ,

where  $\nu$  is unit vector of external normal to the border  $\partial G$ . In particular, at $m=3$ , it should be

${\frac{\partial }{\partial z}\Phi (t,r)|}_{z=0}=0$ .

Algorithms of solving of the initial-boundary-value problems

For solving the considered initial-boundary-value problems one of the following three algorithms is offered to the user at choice in the considered software:⁵

1. Using the classical explicit different scheme;
2. Using the classical implicit different scheme;
3. Using the method of decomposition of the operator.

The problems of optimum choosing of the algorithm parameters on which depend the accuracy, the time and the possibility of practical realization of the equation solution are realized in these algorithms.^5,8

The method of decomposition of the operator is the original algorithm of solving of the diffusion equation use in multidimensional problems. In this algorithm solution of multidimensional diffusion equation is represented in the form of a linear combination of solutions of some one-dimensional diffusion equations.^5,7,8 Thus, the multidimensional problem is reduced to the one-dimensional one. As the experience shows, the algorithm of decomposition of the operator can be realized by computer much more quickly than the classical explicit and implicit different schemes.

Description of river banks by splines

The problem of the analytical description of plane or spatial region for which the diffusion equations and the boundary conditions are investigated is solved.^5,9 This region represents the part of the channel of the river filled with water in the section for which we are modeling water pollution. The interpolation by splines is used for the analytical description of the plane curve, represented in the form of a sequence of distinct points with given Cartesian coordinates. Such a sequence, in particular, could be one of the river bank lines. The construction algorithms for splines of two types of the simplest explicit form which ensure continuous dependence of the tangent vector of the curve on its parameter are realized in the package. In particular, the river bottom is described by polynomial splines, and the river banks are described by trigonometrical splines.^9,10

Mathematical Models Describing the Pollutants Transport in Rivers

The problem of optimum choice of the steps of digitization of the algorithms of all realized different schemes is solved in the considered software.^5,10,11 The criteria of optimality are the requirements of reduction of the time and the errors of calculation as much as possible. The algorithms are proper for the functions describing the polluting substances transport in water having the derivatives up to the fourth order inclusive.

The total number of nodal points  $N$  in the difference schemes determines the time necessary for realization of the algorithm; the accuracy of the obtained result depends on it. Naturally, there arises the question: how to select the numbers $n_1$ , ..., $n_m$  at the given value of  $N$  so that the algorithm should be somewhat optimum.

This problem is solved in the following way:¹¹ at the given values of total nodal points  $N$ , there are determined such values of spatial steps of the grid along coordinate axes  $h_1$ , ...,  $h_m$  for which the upper bound of the module of the residual takes on the minimum value:

$\begin{array}{l}{\displaystyle \sum _{k=1}^m{\varrho }_k}{h}_k^2\to \min ,\\ {\displaystyle \prod _{k=1}^m{h}_k^2}=H=\text{const},\end{array}$

Where  ${\varrho }_k$  is the error of approximation of the equation along the  $k$ th coordinate ( $k=1,\mathrm{...},m$ ).

The solution of this optimization problem is determined by the formula:

$h_k^2=\frac{1}{{\varrho }_k}{\left(H{\displaystyle \prod _{k=1}^m{\varrho }_k}\right)}^{1/m}(k=1,\mathrm{...},m)$

The problem of optimum choice of the value of the parameter  $\tau$  at given spatial steps of the grid is solved in the following way: at the given value of upper bound of the module of the residual of the equation which is equal toε, there are determined such values of  $\tau$  and  $N$  for which the time of calculation takes on the minimum value:

$N^k/\tau \to \min$ ,

$c_{\tau }{\tau }^2+P{N}^{-2/m}=\epsilon =const$ .

Where  ${\varrho }_{\tau }$  is the error of approximation of the equation by the time coordinate;  $k$  is the constant dependent on the algorithm of solution of the sets of linear equations, usually  $k\approx 1$ ;  $\epsilon$  is the given value of the upper bound of the module of the equation residual.

The solution of this optimization problem is determined by the formula:

$N={\left(\frac{P\left(k+1/m\right)}{k\epsilon }\right)}^{m/2}$ ; ${\tau }^2=\frac{1}{c_{\tau }}\frac{\epsilon }{1+km}$ ;

$c_1{h}_1^2=\mathrm{...}=c_m{h}_m^2=k{c}_{\tau }{\tau }^2=\frac{k\epsilon }{1+km}$ .

The main difficulty for practical realization of the described scheme of optimization is connected with the necessity in estimation of the upper bound of partial derivatives of the function  $\Phi (t,r)$  with respect to the spatial coordinates in terms of which parameters  ${\varrho }_{\tau }$  and  ${\varrho }_k$  are expressed without solving the diffusion equation.

 Estimation of derivatives of the unknown function

For practical realization of the described optimization schemes, it is necessary to estimate somehow the upper bounds of the modules of partial derivatives of the unknown function with respect to independent variables without solving the initial task. It should be taken into account that the derivatives of these functions can be estimated with the accuracy of the common constant multiplier.

One way of solution of this problem consists in the following: for some values of the parameters $h_k$  and  $\tau$  (at solving the diffusion equation), the values of the function  $\Phi$  are determined as the first approximation. Then, using the numerical differentiation operations, the required derivatives are determined and the values of the desired parameters are calculated. Using these values, new values of the function  $\Phi$  are determined and so on until the difference between the neighboring calculated values of the function are less than the given value.¹²

At solving the diffusion equation, it is possible to use the explicit scheme at the first stage. Then the determination of the values of the sought for function with higher accuracy is necessary. In this case, for calculation of the unknown function values, the iteration method similar to the one used in work¹ for calculation of the multidimensional integral by the Monte-Carlo method is applied. It is known that the operations of numerical differentiation are not stable, but this should not prevent the realization of the offered method since it requires only rough estimates of unknown derivatives. An obvious drawback of this method is the necessity in performance of a plenty of additional actions. Though, due to the capabilities of modern personal computers, it is negligible at solving the particular tasks.

The effect of smoothness of the inhomogeneous part of the diffusion equation on the accuracy of the results

The inhomogeneous parts of diffusion equation (1) realized in the package contain the impulse functions, which are linear combinations of Dirac delta-functions. These functions and their derivatives are not bounded. Therefore the difference schemes of the solution of diffusion equations are not correct in the vicinities of pollution source localization points.

For elimination of this drawback, it is necessary to consider a model, in which the delta function  $\delta (x-a)$  is replaced by the bounded numerical function  $D(u,x-a)$ , where  $u$  is the additional parameter. It means that the point sources are replaced by extended ones, the capacity of each of which has the maximum corresponding to the capacity of the source at the point of its location. Such sources can be named quasi-point sources.

The function  $D(u,x)$  should satisfy the following requirements: it reaches the maximum value proportional to  $1/u$  at the point  $x=0$ ; it tends to zero at  $x\to \pm \infty$ ; its integral between the limits $-\infty$  and  $+\infty$  is equal to unit; its plot represents a bell-shaped curve. The width of such curve is naturally defined as the width of the rectangle which height is equal to the ordinate of the peak of the curve, and its area is equal to the area of the figure limited by the curve and the axis of abscissas. Thus, the parameter  $u$  characterizes the width of the plot of the function $D(u,x)$ , i.e. the sizes of the source. At  $u\to 0$ , the function  $D(u,x)$ tends to  $\delta (x)$ .

One of the elementary continuous functional dependences which may be used for setting the function  $D(u,x)$  is the trapezoid dependence:

$D(u,x)=\{\begin{array}{ll}1/u\hfill & \text{at }t<1-s,\hfill \\ (1+s-t)/(2us)\hfill & \text{at }1-s\le t\le 1+s\hfill \\ 0\hfill & \text{at }t>1+s,\hfill \end{array},$

Where  $t=2\left|x\right|/u$ ;  $s$ is any real parameter from the interval  $(0,1)$ . The plot of this function together with the axis of abscissas forms an isosceles trapezium, the top and bottom bases of which have the lengths equal to $u(1-s)$  and  $u(1+s)$ , respectively, and the height of which is equal to  $1/u$ . The less is value s, the less this trapezium differs from a rectangle.

Besides the trapezoid dependence, in the developed software package the following dependences are offered to the user’s choice:

1. Gaussian

$D(u,x)=\frac{1}{u\sqrt{2\pi }}{e}^{-\frac{1}{2}{(x/u)}^2}$  ;

2. Lorentzian

$D(u,x)=\frac{2}{2\pi u(1+{(x/u)}^2)}$  .

In fact, the function  $D(u,x)$  defined by one of the two latter relations plays the part of a smoothing function, which allows us to transform any function from  $L^2$  into an infinitely differentiable function by means of the operation of convolution.^10,11

If the pollution sources do not operate continuously, but they operate for a limited interval of time  $[0,T]$ , then the capacity of each source should be smoothed not only by spatial coordinates, but also by time. It means that this function should be proportional to $A(v,t-T)$ , where  $v$  is the additional parameter, and  $A(v,t)$  is the function satisfying the following conditions: at  $t\to -\infty$  it tends to unit, and at  $t\to +\infty$  to zero. Its plot represents a quasi-stepped curve, which steepness is maximum at  $t=0$ , and this maximum value of the steepness is proportional to  $1/v$ . The parameter  $v$  characterizes the time of diminution of the function of discharge. At  $v\to 0$ , the function  $A(v,x)$  tends to  $1-\vartheta (x)$ , where  $\vartheta (x)$  is Heaviside’s stepped function. In the developed software package, it is possible to choose an explicit form of the function  $A(v,x)$ from the following types:

$A(v,t)={\displaystyle {\int }_{t/v}^{\infty }\frac{1}{\sqrt{2\pi }}}{e}^{-{\tau }^2/2}d\tau$ ;

$A(v,t)=\frac{1}{2}-\frac{1}{\pi }\arctan (t/v)$ .

As the computation results show, the smoother are the functions of discharge, the more accurate are the results of solution of the equations in the vicinities of the points of discharge. For given types of functions  $D(u,x)$ and  $A(v,t)$ , the accuracy of the results increases as the values of parameters  $u$  and  $v$  increase.

Input and editing of the initial data

The input of information necessary for operation of the package is performed in the form of separate constants, tables or the choice of the appropriate line from the pull-down. Data for editing are chosen by means of teams of the main menu (Figure 1).

Figure 1 Data for editing are chosen by means of teams of the main menu.

The entered data which form sequences of elements of the same type are represented in the form of tables (Figure 2). When there are several editing windows with tables on the screen, different windows correspond to different sections of the river.

Figure 2 The entered data which form sequences of elements of the same type are represented.

Realization of computation and representation of the results

The results of calculations are displayed in the form of text messages, tables, graphs, etc. For example, the formats of display of computation results as diagrams and tables are shown in Figure 3 & 4, respectively.

Figure 3 The results of calculations are displayed in the form of text messages.

Figure 4 the formats of display of computation results.

The represented text and diagrams may be printed or written in a text- or graphic files, the name of which is determined by default or indicated by the user.

 The software for identification of river water excessive pollution sources located between two controlled cross-sections of the river has been developed by the authors on the basis of mathematical models of pollutant transport in the rivers (described above) and statistical hypotheses testing methods. Depending on the available a priori information, different methods of statistical hypotheses testing can be chosen in the software for identification. With increasing of a priori information, application of more compli­ca­ted methods becomes possible, thus providing their higher assurance. The fol­lo­wing statistical tests are realized in this package:

1. The Euclidean distance method;
2. The Makhalanobis distance method;
3. The unconstrained Bayesian method with the stepwise and arbitrary loss functions;
4. The Constrained Bayesian Method (CBM);
5. The quasi-optimal test based on CBM.

 The peculiarities of CBM methods are described in.¹³⁼¹⁸ Therefore we will not spend space here for their description. Other methods are widely known and described in many scientific works. Let’s show only some characteristics of the package. In particular, in Figure 5 are shown the modes realized in the package, the input of the initial data and the calcu­la­tion regimes. In Figure 6 is shown the graphical view of package output information where the 3-rd point in the bed of the river corresponds to the pollution source guilty of the excessive pollution (Figure 6).

Figure 5 some characteristics of the package.

Figure 6 the graphical view of package output information where the 3^rd point in the bed of the river corresponds to the pollution source guilty of the excessive pollution.

The results of detailed experimental research of algorithms and appropriate programs included in computer packages described ­above are given in.⁵ There are considered the results of research of algorithms of calculation of concentration of polluting substances in the rivers with the help of diffusion equations and also sensitivity of models of different dimensions to the geometrical sizes of a section of the river and dependence of quality of identification ­of emergency pollution sources on the level of the noise, deforming results of measurement. They show truth and reliability both the created algorithms and programs of their realization.^19-23

From the results of calculations, the following practical recommendation concerning the use of the de­ve­­lo­ped models may be given. If geometric sizes, locations of pollution sources, hydrologic characteristics and pollution conditions of the river are such that full mixing of water takes place upstream of the lower controlled section, then it is enough to use the one-dimen­si­onal river models, which are considerably faster and require less computer me­mo­ry than the models of greater dimensionality. Otherwise it is necessary to use the models of greater dimensionality. When choosing a working model for a concrete section of a certain river, it is necessary to per­form preliminary studies with due regard for the above factors. The proper choice of the model and its parameters ensures the qualitative modeling and identification of excessive discharge sources.

The computer packages of realization of mathematical models of propagation of polluting substances in rivers and for identification of exce­ssi­ve pollution sources of the river have been created. The packages comprise original models, methods and algorithms developed by the authors. The software packages are realized for IBM-compatible personal computers according to the generally accepted standard for similar production all over the world. The consumer can use them as a modern, convenient, simple and reliable tool for resolution the problems he deals with in the considered area. Versatile experimental investigation of the developed software packages and the algorithms realized in them, have confirmed their high computing, operational and service qualities. They are good tools for the experts of different professions for qualified solving many practical problems.

None.

The author declares no conflict of interest.

1. Primak AV, Kafarov VV, Kachiashvili KJ. System analysis of air and water quality control. Ukraine: Naukova Dumka: Kiev; 1991.
2. Karaushev AV. River hydraulic. Canada: Leningrad: hydrometeoizdat; 1969.
3. Nikanorov AM, Trunov NM. Intrareservoir processes and quality control of natural water. Hydrometeoizdat: St.-Petersburg, Canada: 1999.
4. Gallakher L, Khobbs DD. Dissemination of pollutions in estuary: mathematical models of the control of water pollution. Moscow: Mir; 1981. p. 228–243.
5. Kachiashvili KJ, Melikdzhanian DI. Advanced modeling and computer technologies for fluvial water quali­ty research and control. New York: Nova Science Publishers; 2012.
6. Gordeziani DG, Gordeziani EG, Kachiashvili KJ. About some mathematical models of propagation of solutions in rivers. IX intergovernmental conference Problems of an ecology and exploitation of objects of power engineering. Sevastopol; 1999. p. 77–79.
7. Kachiashvili KJ, Gordeziani, DG, Melikdzhanian DI. Mathematical models of pollutants transport with allowance for many affecting pollution sources. Urban drainage modeling symposium, Florida: ASCE Library; 2001. p. 692–702.
8. Gordeziani DG, Samarskii AA. Some problems of the thermo-elasticity of plates and shells, and the method of summary approximation. Moscow: Complex analysis and its applications; 1978. p. 173–186.
9. Kachiashvili KJ, Melikdzhanian DI. Analytical description of the bank line of the river for simplification and improvement of the process of calculation of polluting substances concentration. In: Reports of enlarged session of the workshop of I. Vekua Institute of Applied Mathematics. 2002;17(3):101–109.
10. Kachiashvili KJ, Gordeziani DG, Lazarov RG, et al. Modeling and simulation of pollutants transport in rivers. Applied Mathematical Modelling. 2007;31(7):1371–1396.
11. Kachiashvili KJ, Melikdzhanian DI. Parameter optimization algorithms of difference calculation schemes for improving the solution accuracy of diffusion equations describing the pollutants transport in rivers. Applied Mathematics and Computation. 2006;183(2):787–803.
12. Kachiashvili KJ, Melikdzhanian DI, Prangishvili AI. Computing algorithms for solutions of problems in applied mathematics and their standard program realization: part 1-deterministic mathematics. New York, USA: Nova Science Publishers; 2015.
13. Kachiashvili KJ, Mueed A. Conditional bayesian task of testing many hypotheses. Statistics. 2013;47(2):274–293.
14. Kachiashvili KJ. Investigation and computation of unconditional and conditional bayesian problems of hypothesis testing. ARPN Journal of Systems and Software. 2011;1(2):47–59.
15. Kachiashvili GK, Kachiashvili KJ, Mueed A. Specific features of regions of acceptance of hypotheses in conditional bayesian problems of statistical hypotheses testing. Sankhya A: The Indian Journal of Statistics. 2012;74(1):112–125.
16. Kachiashvili KJ, Hashmi MA, Mueed A. Sensitivity analysis of classical and conditional bayesian problems of many hypotheses testing. Communications in Statistics -Theory and Methods. 2012;41(4):591–605.
17. Kachiashvili KJ, Hashmi MA. About using sequential analysis approach for testing many hypotheses. Bulletin of the Georgian Academy of Sciences. 2010;4(2):20–25.
18. Kachiashvili KJ, Hashmi MA, Mueed A. Quasi-optimal Bayesian procedures of many hypotheses testing. Journal of Applied Statistics. 2013;40(1):103–122.
19. Kachiashvili KJ. Comparison of some methods of testing statistical hypotheses: part i. parallel methods. International Journal of Statistics in Medical Research. 2014;3(2):174–189.
20. Kachiashvili KJ. Comparison of some methods of testing statistical hypotheses: part ii. Sequential methods. International Journal of Statistics in Medical Research. 2014;3:189–197.
21. Kachiashvili KJ. The methods of sequential analysis of bayesian type for the multiple testing problem. Sequential Analysis. 2014;33(1):23–38.
22. Kachiashvili KJ. Investigation of the method of sequential analysis of Bayesian type. Journal of Advances in Mathematics. 2014;7(4):1367–1380.
23. Kachiashvili KJ. Constrained bayesian method for testing multiple hypotheses in sequential experiments. Sequential Analysis: Design Methods and Applications. 2015;34(2):171–186.