Software for Pollutants Transport in Rivers and for Identification of Excessive Pollution Sources
Transport and diffusion equation
The diffusion of polluting substances in rivers is mostly described by the threedimensional equation of turbulent diffusion of nonconservative substances [14]:
$\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 timeaveraged concentration of nonconservative 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 timeaveraged speed of the river flow;
$\zeta (r)$
is the coefficient of nonconservatism of pollutants,
$f(t,r)$
is the power of polluting sources.
For the engineering practice very often the twodimensional or onedimensional models are considered instead of the threedimensional one [2,5]. If the nonuniformity of distribution of concentrations of pollutants on the depth of the watercourse is not taken into account, then it is possible to obtain the twodimensional turbulent diffusion equation. The equation of onedimensional 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 [1].
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 onedimensional model).
The additional conditions are specified in the form of
$\Phi (0,r)={S}_{0}\text{\hspace{0.17em}},\text{\hspace{1em}\hspace{1em}}\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 nonclassical. The classical conditions (condition of complete mixing) look like
$\frac{\partial}{\partial x}\Phi (t,r){}_{x=L}=0$
; (condition of full mixing)
nonclassical conditions –
$\Phi (t,r){}_{x=L}=q\xb7\Phi (t,r){}_{x=Ll}$
. (not local boundary condition)
Here
$q$
is the coefficient of selfcleaning of the river on the considered section,
$l$
is the length of the section [5,6].
At using of two or threedimensional models, the following boundary conditions on the other part of the line or on the surface
$\partial G$
(Neumann condition), are set also:
$(\nu \xb7\text{}\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 initialboundaryvalue problems
For solving the considered initialboundaryvalue problems one of the following three algorithms is offered to the user at choice in the considered software [5]:
 Using the classical explicit different scheme;
 Using the classical implicit different scheme;
 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 onedimensional diffusion equations [5,7,8]. Thus, the multidimensional problem is reduced to the onedimensional 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 [11]: 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 \mathrm{min}\text{\hspace{0.17em}},\\ {\displaystyle \prod _{k=1}^{m}{h}_{k}^{2}}=H=\text{const}\text{\hspace{0.17em}},\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}\text{\hspace{1em}\hspace{1em}}(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 \mathrm{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 [12].
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 [1] for calculation of the multidimensional integral by the MonteCarlo 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 deltafunctions. 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 (xa)$
is replaced by the bounded numerical function
$D(u,xa)$
, 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 quasipoint 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 bellshaped 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}t1s,\hfill \\ (1+st)/(2us)\hfill & \text{at}1s\le t\le 1+s\hfill \\ 0\hfill & \text{at}t1+s,\hfill \end{array},$
Where
$t=2\leftx\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(1s)$
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:
 Gaussian
$D(u,x)=\frac{1}{u\sqrt{2\pi}}{e}^{\frac{1}{2}{(x/u)}^{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,tT)$
, 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 quasistepped 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}\mathrm{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.