2020 Математика и механика № 64
УДК 517.3 MSC 93А30
DOI 10.17223/19988621/64/4
V.V. Churuksaeva, A.V. Starchenko
NUMERICAL MODELLING OF POLLUTION TRANSPORT IN TOM RIVER1
This work constructs a mathematical model and a computational method to get extensive data about the structure of a river stream essential for predicting the behavior of a river. The model proposed is based on depth-averaged Reynolds-averaged Navier-Stokes equations. The solver is based on the finite volume method on a staggered structured grid and a high-order Monotonic Upwind Scheme for Conservation Laws for convective fluxes. The solution is obtained with a Semi-Implicit Method for Pressure Linked Equations iterative algorithm based on coupled correction of the depth and velocity fields at each time step. The principal innovation of the algorithm proposed is accounting for the variability of the water depth in the source term in the momentum equations. The results show that the approach proposed accurately predicts the flow field and concentration field and demonstrate the significant role of the flow turbulence in transport of pollution in the river stream.
Keywords: open channel flow turbulence, turbulence simulation and modelling, RANS models, shallow flows, streams and rivers, water quality.
Shallow water equations (SWE) are a widely used approach to modelling geophysical flows that combines acceptable computational cost and precision of the results obtained. It has wide application to modelling problems from fast-developing flows in dam breaks and tsunami waves to relatively slow flows in river estuaries and transport of sediments in lowland rivers.
SWE [1] are obtained by integrating RANS equations by depth and correctly describe a flow in conditions of the depth being much less than the horizontal dimensions of an investigated area. In averaging, the distribution of pressure is hydrostatic, and flow characteristics are assumed to vary little with depth. In this approach, gravity, bottom and free surface shear stresses, and the Coriolis force are the main forces that drive the flow. Because of the assumptions above, SWE are widely used for modelling the atmosphere and the flow in natural basins.
SWE are used to model flow in lakes and seas when the depth-dependent processes do not play a significant role. For example a series of works [2-4] illustrate the application of SWE to modelling tidal flow and tsunami waves in tidal regions of an ocean. Other works [5-7] show the application of this approach to modelling flows with significant wind shear on the surface of the Sea of Azov and the Caspian Sea.
Applications of SWE to modelling flow in rivers and estuaries are extensive [8-10, 3, 11] and are not limited to two-dimensional formulations. 2D SWE are used to compute the flow in river sections that are very long (from dozens of meters to several kilometers) and have an irregular shape that significantly affects the flow [8, 11]. The geometry of the river is especially important in solving such problems as the transport of pollution in the flow, bed deformation [12, 13], and sediment transport [4]. Nevertheless, applications of
1 The reported study was funded by Russian Foundation for Basic Research under research project N 18-3100386.
2D equations to problems of bed deformation and sediment transport are not common due to the very large time scales of these processes. Two of the approaches to modelling processes in a large spatial or time scale are 1D formulations [14-16] and combinations of 1D and 2D models [17, 18] in which 1D computation of the whole river is adjoined with precise 2D computations of the areas of particular interest. 1D solvers can provide accurate information in terms of total sediment load, though they fail at predicting the bed erosion and calculating the final shape of the cross section.
One of the methods to construct the 1D river stream model from 2D depth-averaged equations is integrating them by river width, which gives equations an extra integral parameter—an area of water section. Another approach is omitting the terms that describe the variation of the parameters with one of the spatial coordinates. 1D models are also used when variation of the flow parameters from the river width is negligible compared to longitudinal variation. Very fast-developing flows (such as dam-break problems) or flows in long industrial channels with smooth walls and a flat bottom are examples of flows with such parameters [19].
Solving the fully 3D flow field with specific treatment of the free surface and bottom boundary condition is limited to short river sections adjacent to hydraulic facilities and small-scale flows in laboratory statements [9, 20-22].
In this work, 2D SWE are considered and used as a preferable approach that combines acceptable computational cost for computing long river sections (dozens of kilometers) with the precision of results obtained that is essential for solving such problems as evaluating the anthropogenic changes in the river bed, modelling ice movement and local flooding in spring, and detailing the distribution of pollutants discharged into the river with wastewater.
In many cases, the turbulence of the flow is not accounted for at all (the second derivatives in momentum equations are omitted) [10] or the turbulent viscosity is set as a constant. This approach gives accurate results for cases without recirculation zones or where turbulence is considered mainly to account for energy losses. An example of a flow with such features is coastal waves.
But when turbulent mixing is significant, varying turbulence characteristics should be defined at each point of the flow. To close the SWE, a parabolic eddy viscosity model, a modified mixing length model, and a series of modifications of classical two-parameter difference models such as k-e [23, 24], k-ra [25] have been developed.
The difference in velocity fields obtained for two natural rivers by a depth-averaged solver with standard [26], non-equilibrium [27], and RNG [28] versions of the k-e model is not significant [29]. Turbulent viscosity defined with non-equation models significantly differs from that obtained with the standard, non-equilibrium, and RNG k-e models, which give very close values. In [24, 30] it is shown that the difference between standard models and the Chu and Barbarutsi modification [31] becomes obvious only for flows mainly driven by bottom friction. Simple turbulence models such as the parabolic eddy viscosity model and the mixing length model that are widely used with depth-averaged equations [8, 32] due to their simplicity tend to under-predict the turbulent eddy viscosity values for flows where turbulence is mostly two-dimensional, as it is in river flows [11].
The aim of the work presented in this article is improving the mathematical model and the numerical method constructed in [33] for turbulent river flow computations by adding the modified Streeter-Phelps [34] model of self-cleaning mechanism and illustrating the capabilities of the model computing several test-cases of pollutant transport in Tom river.
Problem Statement
Steady turbulent flow of the viscous incompressible liquid in an open river bed is considered. The river bed has a complex geometry. Assuming that pressure distribution is hydrostatic and water depth is much less than the horizontal size of the area investigated, this flow can be described with 2D SWE:
d(hU) + d(hv) = 0
dx dy
d(hU2) + d(huv) = _gh d(z, + h) +1 d(hXxx) +1dihXxy) + (txz), -(tc)b _ - ,
dx dy dx p dx P dy p x'
. , ^ - Tyx) 1 d (hTyy) (Tyz )s _ ('yz>b „
- + —-- = _gh——-- +--— +-->2— + ^-L--F
d(huv ) + d(hv 2) =_ gh d( zb + h) +1 d (h V ) + 1 d (h Tyy ) + (t yz )s _ (t yz )b _
dx dy dy p dx p dy p y
where h(x,y) is the water depth; u (x,y),v (x,y) are the depth-averaged horizontal velocities; zb (x, y) is the bed elevation; p is the density; g = 9.81m / s2 is the gravity acceleration; Txx, t , xyx, t are the depth-averaged components of the viscous stresses and Reynolds stresses tensor; (x^ )s, (xyz )s, (xxz )b, (xyz )b are corresponding wind stress and bottom friction; and Fx, Fy are the depth-averaged Coriolis force components. The components of the Coriolis force are defined as follows:
_4n _ 4n _
Fx =-hv sin 9, F =-hu sin 9.
time time
Here 9 is the geographical latitude, time is the length of a day in seconds.
In the cases studied, wind stresses (xxz )s, (xyz )s are omitted because their influence
is not relevant compared to the effects produced by the slope and bottom friction terms.
Turbulence Modelling
A Boussinesq hypothesis is used to connect Reynolds stresses with components of the strain velocity tensor by defining effective viscosity (v + vt).
_ , — du dv ^ _ „ , _ .du 2 — _ „ , — ,dv 2 —
=p(v+vt J; Txx = 2p(v+vt ^ - 3 k; = 2p(v+vt ^ - ?k;
(Txz )b = C/ Hu; (yz )b = C/ wlv; C/ = -gf/?.
Here k is the depth-averaged turbulence kinetic energy; v is the kinematic viscosity; vt is the eddy viscosity; cf is the bottom friction coefficient; and n is the Manning
coefficient.
In order to account for production, transport, and dissipation of turbulence in the river flow, a depth-averaged high Reynolds (vt » v) k - e model is used in this work. This model was constructed by Rastogi and Rodi [26] from the original k - e model proposed by Launder and Spalding [35] and was the first differential turbulence model modified for depth-averaged equations.
The equations of the model used are:
vt = S — > ^ e
^^ 4Ih^f l+l(*^f ] + « + P"
dx dy dx y ak dx J dy y ak dy J
d(hue) d(hvs) = d ( vt del d ( vt del ( s„ „ e2 ^
dx
dy
. h^— l+—I h—— 1 +
dx y CTe Sr J cy y CTe dy J
C kPh + Pev - C2 Y
where
Ph =v
2
dv l (du dv
' dy J y dy dx
CTe = 1.3; CTk = 1.0; c1 = 1.44; c2 = 1.92; c„ = 0.09;
Pkv = ck , ; Pev = ce , 2 ; ck = h h2
1
fif
; ce = 3.6-
3/4
f
Here k is the depth-averaged turbulence kinetic energy; e is the energy dissipation; vt is the eddy viscosity; and n is the Manning coefficient.
The constants of the model are assumed to have the same values as in [23] where considerations about their values are also given.
2
k
Boundary Conditions
Inflow boundary conditions are obtained from free surface level and relief data. At the inlet boundary, longitudinal velocity is set to a constant obtained from empirical data for the total discharge of the flow. Because water depth is equal to zero on the river boundary, no wall functions are used and no-slip and no-flow conditions for velocity components are used at wet-dry boundaries. Turbulence kinetic energy and its dissipation fluxes are set to zero at the boundary faces. At the outflow boundary, simple gradient conditions are set for both velocity components and turbulence parameters.
Pollutant Transport Modelling
Water in natural basins has important self-cleaning mechanisms. Wastewater discharged into rivers is diluted by pure river water and partially precipitates. Organic matter introduced into the river with wastewater oxidises by dissolved oxygen driven by microorganisms and algae, and it decomposes from the sun's radiation. The intensity of self-cleaning processes depends on the water level of a river, flow velocity, and the intensity of mixing, which is mainly defined by turbulence. The model proposed contains the following convection-diffusion equation in order to model transport of pollutants whose velocity is equal to the velocity of the flow and whose concentration is relatively small.
dhL + _d_ dt dx
hu L - h I —+ - t
Sc Sct J dx
dL
d_
dy
hv L - h
L Sc Sct J dy
= -k1 Lh - k3 Lh.
Here L (x, y, t) is the depth-averaged concentration of organic matter; k1 is the deoxy-genation rate; and k3 is the rate of organic matter sedimentation.
The problem considered assumes that organic matter is present in the river and that its concentration and the water temperature change slightly with water depth.
The modified dissolved oxygen sag equation [41] is also introduced to the model to define biochemical oxygen demand (BOD) as a criterion of water quality:
dhD d -+ —
dt dx
hUD - h| — + -Vt-
Sc Sct
dD dx
d_
dy
hvD - h | — + -Vt-
Sc Sct
dD
dy.
= k1 Lh - k2 Dh,
where k2 is the reaeration rate and D(x, y, t) is the depth-averaged oxygen deficit, which is defined as the difference between the dissolved oxygen concentration at saturation and the actual dissolved oxygen concentration. In this research, the constants are
set to: k1 =
0.3
k2 =
1.0
; k3 = 0 [36].
In this approach, oxygen deficit is a criterion of both the intensity of self-cleaning processes and water quality.
Both L ( x, y, t )and D( x, y, t ) fluxes are set to zero at the solid boundaries and are set to constants obtained from empirical data at the inflow boundary. At the outflow boundary, simple gradient conditions are used for both variables.
This description of a self-cleaning process is a modification of the classical Streeter-Phelps model. In this model, organic pollution of a river is estimated by BOD and the countervailing influence of atmospheric reaeration. This model could be successfully applied for a timely evaluation of the water quality for a season of the year in several hundred kilometers of a river.
Numerical method
The solver uses a staggered structured mesh to discretise the spatial domain (Fig. 1). Velocity components are defined in the nodes (rhombs) placed in the midpoints of the edges of the initial mesh. All scalar characteristics of the flow are defined in the nodes (circles) of the initial mesh.
■N
n
'ww W 'w
h,
zb u » #
P e
s
S
E
Fig. 1. Control volume and mesh stencil
Discretisation of Convective and Diffusive Terms
Convective fluxes in the momentum equations are approximated with a MUSCL scheme [37] that is up to third-order accurate in regions where function is monotonic.
Values are interpolated from the centers of the mesh cells to the midpoints of the mesh edges with linear function as follows:
— Ax +
®« =°P CTe , Ue > 0;
■s: ^ Ax _ ®e = ®E _yCTe , "e ^ 0;
5 5
where = —^(0e), =— ) are the slopes that are limited by the
Ax Ax v !
function ^(9e ). ^(0) is defined as ^(0) = max
0,mini 20,2+0,2
which gives a
numerical scheme that satisfies conditions of Harten's theorem [38]. In this formulation,
g
0e =]f (§- * 0 ); S« =® e _® p .
An upwind scheme is used to approximate convective terms of the turbulence model equations.
Discretisation of the Source Terms
d( zh + h) d( zh + h) Bed slope source terms _gh-and _ gh- in the momentum
dx dy
equations and turbulence generation due to gradients of horizontal velocities in the equations of the k - e model are discretised with the central difference scheme
hE _hP + z-E _zh
. d(Zh + h)
gh—z—
dx
| gh« —E-p-h-E-^, hE andhP > ew > O(intheriver),
Ax
[0, hE or hP < ew (at thesolid boundary nodes)
because of the construction of the mesh and to obtain second-order spatial accuracy.
The source terms in the equations of the k - e model are represented as S = S - S2®, S > 0, S2 > 0 , where ® is a scalar characteristic (k or e ). This formulation ensures that the resulting value of ® is non-negative.
Computational Method
Equations of the numerical model are solved with an iterative algorithm of coupled correction of the velocity and depth fields that is based on the procedure by Patankar and Spalding for the Reynolds equations [39]. The algorithm used here has been discussed in detail in [33]. The principal innovation of the algorithm is accounting for the variability of the water depth in the source term in the momentum equations. It allows avoiding spurious oscillations and obtaining the correct numerical solution for the fields of velocity and depth. The solver is based on the finite volume method and fictitious areas method and excludes dry cells from computation.
Results and Discussion
The great complexity of modelling flow in a river originates from the complexity of the hydrodynamic model itself, the spatial complexity of the computed domain boundary, the presence of open boundaries, and flow turbulence that affects all processes in a natural basin. Another difficulty is connected with processing the experimental data
about the river bed and general hydrological data about the river, which are essential for initialising the model.
Dommel River Flow Modelling
The mathematical model and the method proposed have been applied to modelling steady turbulent flow in a small shallow river with sharply curved banks . The 250 m section of the Dommel River situated near the border of Belgium and the Netherlands was considered. This object was chosen because the geometry of its bed suggests the formation of large two-dimensional turbulent structures. Another reason is that accurate measurements of its depth, velocity, and bathymetry [40] exist and enable evaluating the results computed.
The bed of the Dommel River is approximately 6 m wide. The flow velocity is set to 0.85 m/s at the inlet boundary, and the water depth is set to 0.3 m at the outlet boundary. The Manning coefficient is set to 0.02 [41]. A structured uniform mesh with 887x401 nodes is used for computations. The convergence of the iterative process is controlled with the value of water discharge at the outlet of the domain. Figure 2 shows the geometry of the studied section of the Dommel River with the traces where measurements of the depth, velocity, and bathymetry were made.
-40
0 -
1.2
0.9
— 0.6
— 0.3
0
220 200 180 160 140 120 100 80 60 40 20 Fig. 2. Geometry of the studied section of the Dommel River (color shows the depth of the river)
Detailed comparison of the numerical and measured velocity fields was made in the cross sections shown in Fig. 2.
As shown in Figure 3, computed and measured values of the depth-averaged flow velocity are in reasonably good agreement within slightly uncertain input data. In the straight parts of the section investigated, velocity profiles are similar to profiles obtained for an open channel accounting for bottom friction and are in good agreement with the data measured. Velocity profiles for cross sections 21 and 34 that are situated in the bends show that the highest velocity magnitude is observed closer to the outer bank of the river and that it gradually decreases toward the inner bank. Recirculation regions form near the inner bank in both bends. The difference between measured and computed velocity values near the banks is connected with the uncertainty in the data about the width of the river. Therefore the results of computations for the S-shaped shallow river show good agreement with measured values and with general concepts.
u, m/s
u, m/s
Fig. 3. Velocity magnitude in several cross sections of the Dommel River (y is distance from the left bank of the river)
Tom River Flow Modelling
The main object of the numerical investigations presented in this paper is the flow in the Tom River (Tomsk Region, Russia).
y, km 1
Basaridaika* a
8 -^Hj^^ s-^-Lwi
Kirgizka
4 i
2
0 2 4 6 8 10 12 14 16 18 20 x, km Fig. 4. Studied section of the Tom River
The Tom River is a representative example of a lowland river; its approximate velocity is 0.33 m/s, and total discharge at the beginning of the investigated section is about 750 m3/s (July). The maximum depth in this section is about 8 m and the mean
-1--.--i-1-1-1-1-1-1—
Basardaika a
-rffc^5^ Ushaika
1 1 »
^^^ Kirgizka
depth is about 2.5 m. The length of the section investigated is 21 km, and the mean width of the river is 800 m. Within the section the Tom River has three tributaries: the Basandaika River (flow rate 2.34 m3/s), the Ushaika River (flow rate 4.35 m3/s), and the Kirgizka River (flow rate 0. 6 m3/s).
In this work, SRTM3 satellite relief data is used to represent the relief of the region. Vertical resolution of the data is 1 m, while the cell is approximately 92x51 m. Local B-spline interpolation [42] is used for data smoothing and defining height in the nodes of the computational mesh.
Critical problems of the river near the city of Tomsk are:
- local flooding in spring connected with ice jams, which occurs annually because of the morphological characteristics of the river [43];
- low water quality near Tomsk because the river suffers from anthropogenic pressure that prevents self-cleaning mechanisms. The city of Tomsk is situated on the bank of the section and organic pollutants are introduced to the river with its drainage waste-water. One of the most considerable organic pollutants introduced into the river in this way is petrochemicals [44]; and
- deformation of the river bed due to anthropogenic and natural factors [43].
Detailed information about the velocity field, depth, and turbulence characteristics is
essential to apply mathematical modelling to solving these problems. Part of the research presented here is applying the numerical model constructed to estimating the water quality in the Tom River.
Transport of pollutants in a river stream is mainly defined by the velocity field (Fig. 5) and the eddy viscosity field (Fig. 6), which defines the intensity of the diffusion.
y, km
8 -
6
4 -
2
0 2 4 6 8 10 12 14 16 18 x, km Fig. 5. Velocity magnitude
y, km 8
6 -
4 -
0 2 4 6 8 10 12 14 16 18 x, km Fig. 6. Eddy viscosity field
The river has a one-branch bed with a small islands in the section studied. The velocity significantly differs from the mean value of 0.33 m/s only in the part between two sharp bends as the bed narrows there. The velocity in the branch between kilometers 14 and 18 is small and pollutants can accumulate in this part of the river. Sharp bends, islands, and inflows situated in the studied area can produce recirculation zones that are the areas of particular interest in this study because they drive the bed deformations and affect the distribution of pollution.
Instant Source Emission Modelling
A case when an organic pollutant has been discharged for 600 s through the inflow boundary in a concentration of 10 mg/l for 600 s has been calculated to simulate a situation when a pollutant has escaped into a small area in a river from an industrial facility failure. Figure 7 illustrates the movement of the polluted water transport and dispersion in a river stream within 30 hours after discharge through the inflow boundary and shows how pollution is transported by the flow and dispersed for 30 hours after the discharge.
y, km
86
42
0 2 4 6 8 10 12 14 16 18 x, km Fig. 7. Transport of the pollution from the instant source
The pollutant is transported by the stream and at the same time its maximal concentration decreases with time because of mixing and oxidizing of the organic substance; pollution is removed from the section computed in approximately 30 hours (Figure 7).
Dispersion of the pollution is maximal at the narrowest section between two bends where large eddy viscosity causes significant turbulence mixing, and convective flux is large due to the highest velocity. Distribution of the pollutant across the river is almost uniform in the initial section. Moving downstream, the distribution substantially expands following the shape of the velocity profile.
The following diagrams show a detailed illustration of the transport of pollution and oxygen deficit variation (for the case with self-cleaning processes considered) along the river with time to estimate the influence of the self-cleaning processes to the water quality in the river.
Maximal concentration computed for the case with the self-cleaning mechanisms considered is lower than concentration computed for the second case. In both cases, concentration of the pollutant significantly decreases (from 1 mg/l to 0.36 mg/l) moving from kilometer 7 to 15 downstream from the initial cross section (Figure 8).
1
0.8
0.6 0.4
0.2
0
0.05 0.04 0.03 0.02 0.01
0 5 10 15 x, km
Fig. 8. Variation of the width-averaged values of concentration and oxygen deficit with time
The initial value of the oxygen deficit is set to zero in order to consider that the concentration of the dissolved oxygen is enough to initiate a self-cleaning process. With oxidizing, the deficit increases and reaches a maximum value of 5 mg/l at t = 8 h when the oxidation process is the most intensive and then decreases as a result of the decrease in the concentration of the organic matter in the water and reaeration.
Permanent Source Emission
The tributaries of the Tom River flow through the urbanised territory and collect drainage water with a high concentration of organic pollutants (phenol and petrochemicals). The case with permanent sources of pollution at the mouths of the Ushaika and Basandaika rivers was computed to investigate distribution of the pollution that comes to the Tom River from its tributaries. In this case concentration of the pollutant discharged is 10 mg/l.
The intensity of the pollution discharge from the Ushaika is more than that from the Basandaika because its flow rate is larger. The distribution of pollutant near the discharges is significantly non-uniform along the cross section of the river because the tributaries' flow rates are much less than the flow rate in the Tom. Pollution is transported mainly along the bank of the river where tributaries inflow. Near the inflows, concentration of the pollutant is maximal. Moving downstream, pollution is dispersed more uniformly due to the significant turbulence mixing in the river stream.
Pollutant concentration The case with self-cleaning mechanisms (solid line). The case without self-cleaning mechanisms (dashed line).
----4 hours
----8 hours
16 hours 24 hours
29 hours .....
after dischargement
Oxygen deficit The case with self-cleaning mechanisms
-1-1-1-1-1-1-1-1-r
4 6 8 10 12 14 16 18 x, km
Fig. 9. Concentration of a pollutant discharged from the inflows
2
Oxygen deficit distribution shown in Figure 10 mainly follows the picture of pollutant concentration - increases in areas with higher concentration due to the intensity of the oxidising process and decreases where there is not much pollutant. Oxygen deficit also high in the recirculation zone after the second bend, indicating a higher concentration of the pollutant there.
To evaluate the influence of the turbulence on pollutant transport, the permanent source emission case was computed without considering turbulence in both the pollutant transport equation and modified dissolved oxygen sag equation. Comparison of the results obtained shows that in the case without the turbulence accounted for, concentration of the pollutant near sources is 15-25% greater because of smaller diffusion in the stream. This quantitative estimate was made in the cross section of the river at x = 12000 m.
Conclusions
Recirculating areas are of particular interest for this research because they tend to accumulate pollution. The mathematical model constructed allows identifying both areas with large convective flux and those with intensive turbulent mixing due to large eddy viscosity. The computed fields of the velocity and depth obtained for the Dommel
River show good agreement with measured values, which confirms that the model constructed correctly represents flow in a curved shallow river. Computed values for the velocity and concentration fields in the Tom River show the significant role of the turbulence in the flow structure and in distributing pollutants. Computed values of the concentration have valuable non-uniform distribution along a cross section of the river, so that a mathematical model with two spatial coordinates should be used to obtain correct and detailed data about the distribution of the pollutant.
With precise bathymetric data, the approach proposed could be used to produce a detailed description of the current structure of the river stream and a quantitative forecast for the distribution of the pollutants in the stream.
Acknowledgements
The authors wish to thank Jean Kollantai of Tomsk State University for her assistance with style.
REFERENCES
1. Kuipers J. and Vreugdenhil C.B. (1973) Calculations of two-dimensional horisontal flow. Delft Hydraulic laboratory, Delft, Report on basic research S 163.
2. Marchuk A.G. (2015) Computing of tsunami heights above the inclined bottom relief within the wave-ray approach. Siberian Journal of Numerical Mathematics. 18(4). pp. 377-388 [in Russian].
3. Sauvaget P., David E., and Soares C. (2000) Modelling tidal currentsonthe coast of Portugal. Coastal Engineering. 40. pp. 393-409.
4. Castro Diaz M.J., Fernandez-Nieto E.D., Ferreiro A.M. (2007) Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods // Computers & Fluids. 37. pp. 299-316.
5. Krukier L.A. (2000) Mathematical modeling of hydrodynamic processes in Azov Sea. Patterns of Oceanographic and Biological Processes in Azov Sea. Apatity: KSC RAS. pp. 129-163 [in Russian].
6. Chikin A.L. (2001) Construction and numerical investigation of the 3D hydrodynamic model of Azov sea. Proceedings of International Conference RDAMM-2001. 6(4). pp. 686-692 [in Russian].
7. Shabas I.N. (2014) Mathematical modeling of the transport of multicomponent contaminant in Azov sea on multiprocessor systems. Izvestiya SFedU. Engineering Sciences. 12(161). pp. 200-210 [in Russian].
8. Yu L., Zhu S.P. (1993) Numerical simulation of discharged waste heat and contaminants into the south estuary of the Yangtze River. Mathematical and Computer Modelling. 18(12). pp. 107-123.
9. Olsen N.R.B., Stokseth S. (1995) Three-dimensional numerical modelling of water flow in a river with large bed roughness. Journal of Hydraulic Research. 33. pp. 571-581.
10. Hou J., Simons F., Mahgoub M., and Hinkelmann R. (2013) A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography. Computer Methods in Applied Mechanics and Engineering. 257. pp. 126-149.
11. Cea L., Puertas J., and Vazquez-Cendon M.E. (2007) Depth averaged modelling of turbulent shallow water flow with wet-dry fronts. Archives of Computational Methods in Engineering. 14(3). pp. 303-341.
12. Zinov'ev A.T., Marusin K.V., Shibkih A.A., Shlychkov V.A., Zatinatskij M.V. (2006) Mathematical modeling of the flow dynamics and bed deformations on the section of Ob' river near the city of Barnaul. Polzunovskij vestnik. 2. pp. 204-209 [in Russian].
13. Zinov'ev A.T., Marusin K.V., Shibkih A.A., Shlychkov V.A., Zatinatskij M.V. (2006) Modeling of bed processes to evaluate results of dredging the river bed. Polzunovskij vestnik. 2. pp. 197-203 [in Russian].
14. Lyubimova T.P., Lepihin A.P., Parshakova YA.N., Tiunov A.I. (2010) Numerical modeling of dilution of highly mineralized brines in turbulent flows. Computational continuum mechanics. 4(3). pp. 68-79 [in Russian].
15. Karepova E.D. (2008) Modeling of the unsteady flow in the lower pool of the Boguchany HES. Computational Technologies. 13(2). pp. 28-38.
16. Belolipetskij V.M., Genova S.N., Petrashkevich V.I. (2001) Numerical modeling of pollutant transport in a river flow. Proceedings of International conference RDAMM-2001. pp. 127133 [in Russian].
17. Finaud-Guyot P., Delenne C., Guinot V., and Llovel C. (2011) 1D-2D coupling for river flow modeling. Comptes Rendus Mecanique. 339. pp. 226-234.
18. Fernandez-Nieto E.D., Marin J., and Monnier J. (2010) Coupling superposed 1D and 2D shallow-water models: Source terms in finite volume schemes. Computers & Fluids. 39. pp. 1070-1082.
19. Bulatov O.V., Elizarova T.G.: (2011) Regularized shallow water equations and an effective method for numerical modeling of the flow in shallow basins. Journal of Computational Mathematics and Mathematical Physics. 1(51). pp. 170-184 [in Russian].
20. Kang S., Lightbody A., Hill C., and Sotiropoulos F. (2011) High-resolution numerical simulation of turbulence in natural waterways. Advances in Water Resources. 34. pp. 98-113.
21. Kang S., Sotiropoulos F. (2012) Numerical modeling of 3D turbulent free surface flow in natural waterways. Advances in Water Resources. 40. pp. 23-36.
22. Sandiv S.K., Sotiropoulos F., Odgaard A.J. (1998) Three-dimensional numerical model for flow through natural rivers. Journal of Hydraulic Engeneering. 124 (1). pp. 13-24.
23. Rodi V. (1984) Models of turbulence in the environment. Methods of Measuring Turbulent Flows. pp. 276378 [in Russian].
24. Babarutsi S., Chu V.H. (1998) Modelling transverse mixing layer in shallow open-channel flows. Journal of Hydraulic Engineering. 7(124). pp. 718-727.
25. Yu L., Righetto A.M. (2001) Depth-averaged k-omega turbulence model and application. Advances in Engineering Software. 32. pp. 375-394.
26. Rastogi A.K., Rodi W. (1978) Predictions of heat and mass transfer on open channels. J. Hydraul. Div. HY. pp. 397-420.
27. Kim S.E., Choudhury D. (1995) A near-wall treatment using wall functions sensitized to pressure gradient. ASME FED. 217, Separated and Complex Flows.
28. Yakhot V., Orzsag S.A., Tangam S., Gatski T.B., and Speziale C.G. (1992) Development of turbulence models for shear flows by a double expansion technique. Phys. Fluids A. 4(7). pp. 1510-1520.
29. Wu W., Wang P., and Chiba N. (2004) Comparison of five depth-averaged 2-D turbulence models for river flows. Archives of Hydro-Engibeering and Environmental Mechanics. 51(2). pp. 183-200.
30. Babarutsi S., Nassiri M., and Chu V.H. (1996) Computation of shallow recirculating flow dominated by friction. Journal of Hydraulic Engineering. 122(7). pp. 367-372.
31. Barbarutsi S., Chu V.H. (1991) A two-length-scale model for quasi two-dimensional turbulent shear flows. Proc 24th congr of IAHR. C. pp. 51-60.
32. River2D Hydrodynamic Model for Fish Habitat [Resource]. River2D: [website]. [2002]. URL: http: //www.river2d.ualberta.ca/.
33. Churuksaeva V.V., Starchenko A.V. (2015) A mathematical model and numerical method for computation of a turbulent river stream. Tomsk State University Journal of Mathematics and Mechanics. 6(38). pp. 100-114 [in Russian].
34. Gromova V.V., Mihailov M.D. (2011) Numerical modeling of self-cleaning processes in the river accounting for additional wastewater treatment. Current problems in mathematics and mechanics: The second state young scientists conference devoted to the 90th anniversary from the day of academician N.N. Yanenko's birthday (October 12-14, 2011). pp. 254-259 [in Russian].
35. Launder B.E., Spalding D.B. (1974) The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering. 2(3). pp. 269-289.
36. Vavilin V.A. (1983) Nelineynye modeli biologicheskoy ochistki i protsessov samoochish-cheniya v rekakh [Nonlinear models of biological wastewater treatment and self-cleaning processes in rivers]. Moscow: Nauka [in Russian].
37. Cada M., Torrilhon M. (2009) Compact third-order limiter functions for finite volume methods. Journal of Computational Physics. 228. pp. 4118-4145.
38. Harten A. (1983) High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics. 49. pp. 357-393.
39. Patankar S.V. (1980) Numerical heat transfer and fluid flow. Washington, DC: Hemisphere Publishing Corporation.
40. de Vriend H.J. and Geldof H.J. (1983) Main flow velocity in short and sharply curved river bends. Department of Civil Engineering Delft University of Technology. Delft. 83.
41. Chow V.T. (1959) Open Channel Hydraulics. New York: McGraw-Hill.
42. Zav'yalov Y.S., Kvasov B.I., Miroshnichenko V.L. (1980) Metody splayn-funktsiy [Methods of spline functions]. Moscow: Nauka [in Russian].
43. Tarasov A.S., Vershinin D.A. (2015) Building a predictive model of ice jam occurrence on the branched site of the tom river with hydrological computer modelling. Tomsk State University Journal. 390. pp. 218-224 [in Russian].
44. L'gotin V.A., Makushin Y.V., Savichev O.G. (2005) Information Resources. TomskgeomonitoringRegional Center. URL: http://www.tgm.ru [in Russian].
Received: December 27, 2019
Vladislava V. CHURUKSAEVA (Candidate of Physics and Mathematics, Faculty of Mechanics and Mathematics, Tomsk State University, Tomsk, Russian Federation). E-mail: chu.vv@mail.ru
Alexander V. STARCHENKO (Doctor of Physics and Mathematics, Professor, Faculty of Mechanics and Mathematics, Tomsk State University, Tomsk, Russian Federation). E-mail: starch@math.tsu.ru
Churuksaeva V.V., Starchenko A.V. (2020) NUMERICAL MODELLING OF POLLUTION TRANSPORT IN TOM RIVER. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 64. pp. 48-62
DOI 10.17223/19988621/64/4