Mathematical Model of Suspended Particles Transport in the Estuary Area, Taking into Account the Aquatic Environment Movement

: A large amount of contaminants enter marine systems with river runoff, so the purpose of the study is to model the transport of suspended particles in the estuary area. To describe hydrodynamic and hydrophysical processes, the mathematical model of the suspended particles transport was used, supplemented by a three-dimensional mathematical model of hydrodynamics, used to calculate the ﬁelds of the aquatic environment movement velocity vector, and equation for calculating the variable density. The approximation of the equations for calculating the velocity ﬁeld by spatial variables is based on the splitting schemes for physical processes with ﬂuid volume of the control areas, which allows for us to consider the complex geometry of the coastline and the bottom. The suspended particles transport model is approximated using splitting schemes for two-dimensional and one-dimensional problems. Numerical experiments were carried out to simulate the aquatic environment movement in the estuary area, the multicomponent suspension deposition, as well as mixing of waters in the mouth, taking into account the different density of the aquatic environment. The used models and methods allow to signiﬁcantly improve the accuracy of modeling suspended particle transport and consider the factors inﬂuencing the studied processes.


Introduction
In recent decades, pollution of the planet's water resources has become increasingly threatening. This situation is mainly related to human activity, which lead to the fact that a significant amount of pollutants enters aquatic ecosystems, and subsequent resuspension, caused by natural or human-made situations contributes to secondary pollution. All this negatively affects the production and destruction processes in aquatic ecosystems. The study of the natural mechanisms of hydrodynamic processes, taking into account the variability of meteorological conditions, makes it possible to choose a strategy to minimize damage from water pollution or prevent pollution of the aquatic ecosystem.
In the Russian Federation, scientific research on the creation and research of complex marine systems mathematical models has more than half a century history. Many scientists actively researched the issue of optimal exploitation of water resources, the development of models for the transport of contaminants in water bodies and the study of assessing their impact on the biological resources of a water body. In the work of Matishov G.G. [1], the water exchange between the Black and Azov seas on the basis of long-term ship observations of the Azov-Black Sea basin is being investigated. In the work of Ilyichev V.G., the evolutionary stable characteristics of the Azov Sea are investigated with variations in averaged model for estuaries in [16]. Chen XJ. proposes the methods for shallow water models, which solve unsteady Reynolds-averaged Navier-Stokes equations. Experiments were carried out with an idealized estuary case which has a large basin and two narrow tributaries. In the work [17] Bradford S.F. describes the model based on Navier-Stokes equations with using of a modified sigma coordinate transformation. This model allows us to simulate free surface flows and is used for simulating fluid-structure interactions.
A study of the current achievements of Russian and foreign scientists in the field of mathematical modeling of aquatic ecology processes has shown that the mathematical models that are part of the software systems used to predict changes in the state of aquatic ecosystems when natural and technogenic phenomena occur in them do not take into account important factors that have a significant impact on the nature of the course of the processes under study. Simplifying hydrodynamic models does not consider turbulent exchange, the Coriolis force, the complex geometry of the coastline and the bottom, wind stresses and friction against the bottom, wind and surge phenomena, evaporation, and river flows. In the coastal areas of shallow water bodies, salt and fresh waters mix and there is a significant difference in depths. The problem of constructing discrete analogs of the developed mathematical models with the property of stability arises when modeling the hydrophysical processes of shallow water bodies. The hydrophysical processes can be described by non-stationary and spatially non-homogeneous mathematical models, including a 3D model of the hydrodynamic processes of a reservoir, models for the transport of salts, heat, and suspension, which can be represented as systems of nonlinear partial differential equations. Such problems can be solved based on the methods for solving diffusionconvection equations. Another urgent task in the construction of mathematical models of hydrodynamics for the prediction of storm events, the transport of pollutants, and other dangerous natural and human-made phenomena is the development of difference schemes under the condition that convective transport prevails over diffusion transport [18][19][20][21]. When modeling hydrodynamic processes in channel systems, large values of grid Péclet numbers arise, therefore, the accuracy of convective transport modeling decreases. It is also necessary to use models that consider changes in the density of the medium due to a significant change in salinity. Such systems are difficult to model, as re-suspension occurs and the structure of the bottom changes dynamically due to sedimentation. Standard difference schemes do not work. In this paper, methods of mathematical modeling are proposed that allow to improve the accuracy of modeling the listed processes.
The purpose of this work is to build a three-dimensional non-stationary hydrodynamics model coupled with the model of multicomponent suspended matters transport, as well as to develop effective methods for the numerical implementation of these problems. In this work, proposed hydrodynamics model which considers evaporation and precipitation as in the continuity equation, and in the motion equations, along with the complex shape of the coastline, the Coriolis force, wind stress and friction on the bottom, which has a complex relief, etc. This model based on 3D Navier-Stokes equations and the continuity equation for an water medium with variable density. The particles diameter, the difference between the particle density and the liquid density and the dynamic medium viscosity consider in the model of multicomponent suspended matters transport. In the numerical implementation of hydrophysics models, which can be reduced to problems of the diffusion-convection type, a developed difference scheme is used. A proposed difference scheme is an optimized scheme based on the Standard Leapfrog and Upwind Leapfrog schemes. The modified difference scheme contains weight coefficients equal to 2/3 and 1/3. When calculating the coefficients, the order of the approximation error was minimized. If the grid Péclet number is sufficiently large, it becomes difficult to construct a difference scheme of a high order of accuracy in the mathematical modeling of hydrophysical processes in reservoirs with complex bathymetry. Such cases arise when modeling water flows in river beds. The method of filling of rectangular cells with a material medium, in particular, with a liquid, is used to improve the smoothness and accuracy of the approximation of solving hydrodynamic problems in a region of complex shape. The above methods help to improve the accuracy of numerical simulation of hydrodynamic processes in a region of complex shape. As an illustrative example of the proposed methods using, the problem of transporting suspended matter from the riverbed to the sea in the estuary areas is solved.

Problem Statement
To predict natural and human-made hazards, system of initial-boundary value problems was constructed, which includes the equations of hydrodynamics, multicomponent suspended matters transport and variable density.

Hydrodynamics Model
The considered mathematical model of the hydrodynamics of the estuary area includes [22]: • equations of motion (Navier-Stokes equations): • continuity equation for variable density: ∂ρ ∂t where V = {u, v, w} is water flow velocity vector for the studied shallow water body where n is the normal vector directed inside the computational domain, ω is the evaporation rate of liquid, and τ x and τ y are components of the tangential stress.
On the free surface, we set the tangential stress as: τ x , τ y = ρ a Cd s |w| w x , w y , where Cd s = 0.0026, where w is a vector of the wind speed relative to the water [m/s]; ρ a is the density of the atmosphere [kg/m 3 ]; Cd s is surface resistance coefficient, which can be set in the range 0.0016-0.0032. For definiteness, we assume that it depends on the wind speed.
The tangential stress vector for the bottom, considering the movement of water, is set as: τ x , τ y = ρCd b |V|{u, v}, Cd b = gn 2 /h 1/3 , where h is the depth of the studied reservoir [m]. The coefficient of roughness n=0.04 is set according to Manning. Consider that n ∈ [0.025, 0.2].

Model of Transport of Multicomponent Suspended Matters
Transport of suspended multicomponent particles is described using the diffusionconvection equation of the following form: The process of particle settling occurs according to the laws of falling bodies in a medium that resists their movement. When settling, the particles first move rapidly, and then the frictional resistance force of the medium and the force of gravity are balanced, and the particles acquire a constant speed and settle evenly.
As parameters that have a significant impact on the rate of sedimentation of suspension, we can distinguish: • particle diameter (the diameter is greater, the particle settling rate is greater); • the difference between the particle density and the liquid density (the difference in densities is greater, the settling rate is higher); • dynamic medium viscosity (the dynamic medium viscosity is lower, the particle settling rate is higher).
The constant settling rate can be determined by the formula (Stokes' law): where d r is the diameter of the deposited r-th fraction particle [m], ρ V,r is the density of the deposited r-th fraction particle [kg/m 3 ], ρ is the density of the medium [kg/m 3 ], η is the dynamic medium viscosity [Pa·s]. It should be noted that the use of Stokes' law is possible only within certain limits. The upper limit is determined by the moment of transition from suspension to colloidal solutions, when the particles of the dispersed phase have a size of 0.1-0.5 η, and also considers the effect of Brownian motion, which does not prevent particle settling.
The upper limit of the use of the Stokes' law is characterized by the numerical indicator of the Reynolds criterion Re ≈ 2. In the event that the resistance of the medium is proportional to the square of the velocity and Re > 2, then the following formula is used to calculate the particle settling velocity: where ζ is the medium hydraulic resistance coefficient. At 2 < Re < 500, the value of the resistance coefficient ζ = 18.5/Re 0.6 , and in the case of 500 < Re < 1500, the resistance coefficient is ζ = 0.44. Almost always, the sedimentation rate in a liquid medium is determined by the numerical value of the Reynolds criterion with a preliminary determination of the value of the Archimedes criterion. Even in coarse suspensions, as a rule, there are a sufficient number of particles for which Re < 2. Thus, they have a low settling rate, which can be determined from Stokes' law.

Equation of Variable Density
The standard formula is used to calculate the density of the aquatic environment: where V r , ρ V,r are the volume fraction and density of r-th fraction of suspension, ρ 0 is the density of fresh water under normal conditions, and R is number of fractions. The density is found from Equation (4). The found solution is substituted into the system of motion Equations (1) and the continuity Equation (2). In the absence of impurity, the continuity equation looks like div(V) = 0. The solution of the system of Equation (1) is found taking into account Equation (2).

Approximation of Hydrodynamics Equations
Introduce a uniform grid for approximation of the 3D hydrodynamics mathematical model (1)- (2): where τ is the time step; h x , h y , h z are the space step; N t is the number of time layers; T is the upper bound in time; N x , N y , N z are the number of nodes in space; l x , l y , l z are the characteristic dimensions of the computational domain. The method of splitting by physical processes is used. In this case, the hydrodynamic model is transformed into three subproblems [23][24][25]. The components of the water flow velocity vector will be calculated on the basis of the first subproblem of the diffusionconvection type selected in the process of splitting. For calculations on the intermediate time layer, we will use the following equations: where V = {u, v, w} is the velocity vector with components in the previous time layer.
In Equation (5), we introduced the notationũ,ṽ,w are the components of the velocity vector in the intermediate time layer To solve the second subproblem (calculation of pressure), the equation is used: whereρ and ρ are distribution of the density of the aquatic environment on the current and previous time layers, respectively. Equality (6) is based on the Poisson equation.
When solving the third subproblem, obtained as a result of splitting into physical processes, the water flow velocity is determined in a new time layer using the formulas: whereû,v,ŵ are the components of the velocity vector at the current time layer. Formula (7) are of explicit type and are easy to implement. The degree of cell "fullness" in the implementation of the calculation algorithms is indicated o i,j,k and is determined on the basis of the pressure of the water column on the bottom of the considered cell (i, j, k). According to the expression from [26], the degree of cell filling is determined as follows: The position of the free surface is determined by the pressure of the liquid column based on Formula (8). If the liquid is not in a state of weightlessness, this method of calculating the level elevation function guarantees the conservation of the amount of substance.
From the total hydrodynamic pressure P(x, y, z, t), two components can be conditionally distinguished the pressure of the liquid column or hydrostatic pressure and the excess pressure over hydrostatic pressure: where ρ 0 gz is the hydrostatic pressure of the undisturbed fluid, p is the excess pressure over hydrostatic pressure. Using (9), three subproblems (5)-(7) will be written in the forms: Note that the term g(ρ 0 /ρ − 1) in the last equation from (10) describes buoyancy or the Archimedes force.
The separation of two components from the pressure is caused by computational expediency. The calculation of problem (11) is computationally expensive because the operator is ill-conditioned. The pressure is calculated in two steps: • first step: the pressure is calculated in the hydrostatic approximation (without taking into account vertical movement and depth stratification) based on two-dimensional subproblem; • second step: the pressure excess over the hydrostatic pressure is found.
Such an approach for calculating the hydrodynamics of shallow reservoirs can significantly reduce computational costs. Note that p, in contrast to P, practically does not change with depth.
Applying the balance method, modified by introducing fill factors, an approximation of the hydrodynamic problem is implemented, as a result of which the calculation of the water flow velocity field is calculated.

Approximation of the Suspended Particles Transport Model
We approximate by an equation of the diffusion-convection type for the threedimensional case: Equation (13) is an analog of Equation (3) for the one-component suspension. (Using the example of Equation (13), an approximation of Equations (3) and (10)-(12) is considered.) A grid (uniform rectangular) ω τ was used for calculations. The time step τ was used for discretization: For Equation (13), the splitting schemes into a two-dimensional and one-dimensional problem are used [27,28]:  (14) and (15), respectively. Solve a model problem of the form (10). When organizing calculations, a grid is used in the form: where h x , h y , τ are the space and time steps, respectively, N x , N y , N t are the upper bounds in time and space, respectively. Denote the characteristic dimensions of the computational domain in the directions of the axes Ox, Oy in the following way: l x , l y .
Set the volume of fluid (VOF) of the areas located near the cell using the coefficients q 0 , q 1 , q 2 , q 3 , q 4 . The filling of the area D 0 is determined by the coefficient q 0 : Denote Ω m m = 0, 4 the parts of the areas D m that are filled. Calculate the coefficients q m based on the works [29]: Let the boundary condition have the form ∂c ∂n = α n c + β n . Discretization the convection and diffusion operators is calculated on the basis of [29]: Approximation of model (14) is based on a difference scheme modified by specifying a combination of Standard Leapfrog and Upwind Leapfrog schemes. In this case, the values of the weight coefficients of the linear combination used in the difference scheme can be calculated on the basis of minimizing the order of the approximation error, as was done in [18,29]. This difference scheme increases the accuracy of the hydrodynamic problems numerical solution for large values of the grid Péclet number (in the range from 2 to 20) [29]. Additionally, the approximation of model (14) considers the function values of filling the cells of the computational domain [18]. Equation (15) is solved by the elimination method [28].
For the function q(x, t), set the sum of the Fourier series (finite, in complex form) according to the formula: where ω = π L , m is the number of the harmonic, C m (t) = 2 L L 0 q(x, t) exp(−jωmx)dx is the complex amplitude of m-th harmonic, j = √ −1. The sum of the components of different frequencies is written on the right side of Formula (17). When constructing a numerical solution, consider a grid (uniform, rectangular), define it as a set of the following form: ω h = x i = ih; i = 0, N; Nh = L . When setting the grid, the following parameters were used: h, N -step and the number of nodes in the spatial coordinate. Discretization of Equation (16) using the introduced grid allows us to pass to an equation of the form: Note that the eigenvalues of the operator (Λq) i can be found by the formula: Consider a space-time grid ω = ω h × ω τ by introducing an additional time grid ω τ = t n = nτ, n = 0, N t , N t τ = T , where τ, N t is the step and number of nodes in time, respectively. To approximate Equation (16) with respect to a time variable, we will use scheme with weight: where C n+σ 1] is the weight of the scheme, χ m = λ m τ. Let us find the value of the weight σ of the scheme, at which the error of the numerical solution is minimal. The value of the error at the n-th time layer, in terms of the exact value C m (t n ) and the approximate value of the function C n m , can be expressed in terms of the following function: The estimate of the error in the numerical solution of the problem (16) based on [18]: It is obvious that the approximation error with respect to the time variable on the i-th layer does not exceed max m |ψ n m |. Let us introduce the notation: Then, taking into account (20) and Parseval's theorem, we obtain: Consider the function: +2(Re z 1,m Re z 2,m + Im z 1,m Im z 2,m )σ + (Re z 2,m ) 2 + (Im z 2,m ) 2 .
To find min σ |s m (χ m , σ)|, we define the extremum points of the function s m (χ m , σ) by the variable σ: From (21) obtain: The above reasoning allows us to determine the optimal value of the scheme weight σ at which the error in the numerical solution of problem (16) is minimal.

Results of Numerical Experiments
Software package in C++ for the numerical solution of problem (1)-(4) was developed. This software considers such physical parameters as: turbulent exchange, complex bathymetry, the influence of wind and friction on the bottom surface of the study area, and the presence of a significant density gradient in the aquatic environment. The software package allows us to calculate three-dimensional water flow velocity fields based on the model (1) and (2), the model of transport of suspended particles (3) and the model of transport of suspended particles, taking into account the movement of the aquatic environment (1)-(4).

Numerical Implementation of the Hydrodynamic Model
In the numerical implementation of the model (1) and (2) based on the developed software package, the following functions associated with the calculation are performed: • water flow velocity fields, while pressure will not be considered; • hydrostatic pressure, which in the numerical implementation of the considered mathematical model can be used as an initial approximation when calculating the values of the hydrodynamic pressure function at the nodes of the computational domain of the considered uniform rectangular grid; • hydrodynamic pressure; • fields of water flow velocity in the three-dimensional case.
Consider modeling the movement of the aquatic environment using a test problem. The uniform rectangular grid 100 × 100 × 40 computational nodes is introduced. The parameters of the computational domain are: 50 × 50 × 4 (length, width and depth are determined). The horizontal steps were 0.5 m, the vertical-0.1 m. The calculations were carried out for a time interval of 1 min with a time step of 0.25 s. The input data for the calculation according to model (1) and (2): pressure is 1.29 Pa, water density is 1000 kg/m 3 , the horizontal and vertical components of the turbulent exchange coefficient are 0.01 m 2 /s and 0.0005 m 2 /s, respectively, and the acceleration of gravity is 9.81 m/s 2 .
Determine the geometry of the computational domain using a depth map (Figure 1). In Figure 1, legend on the right the depth values in meters are shown. Figure 2 shows a numerical experiment based on the developed software package. The color reflects the values of the water flow velocity vector (three-dimensional case). A part of the computational area, the estuary area, has been identified. The calculations were carried out on the basis of the developed software package for various time intervals: 15 s, 30 s, 45 s and 1 min. The considered scenario ( Figure 2) allows us to observe the complex dynamics of the water flow movement process in process of time, which is quite typical for the researched phenomena.
The developed software package allows us to predict the appearance of flooding and drainage areas in the estuary area of the river.

Suspended Particle Transport Modeling
Describe the software implementation of the suspension transport model for the test problem of hydrophysics of the estuary area, in which the process of sedimentation of two fractions with different properties and characteristics is studied. For fraction A, set the deposition rate to 2.4 mm/s. Let the percentage of fraction A in dusty particles be 36%. In this case, we set the sedimentation rate of fraction B to 1.775 mm/s. Let the percentage of fraction B be 64%. Computational domain parameters: length 1 km; width 720 m. In the experiment, the horizontal and vertical steps were 10 and 1 m, respectively. The calculated interval within the considered experiments was 2 h. In the framework of the described experiment, the suspension source is placed at a distance of 5.5 m from the bottom. Determine the average flow velocity in the area with the suspension source, its value is 0.075 m/s. Build graphs of changes in the granulometric composition of the bottom. We will set different initial concentrations of suspended particles ( Figure 3). As part of the experiment, we believe that a horizontal axis passes through the dredging area, directed along the current. On the vertical axis of Figure 3a,b show the depth of the reservoir (m). Let us direct the Oz axis down vertically. The level of bottom sediments (in mm) is plotted along the vertical Oz axis directed vertically upwards (Figure 3c,d). During the experiment, the concentrations of fraction A in water were determined (Figure 3a). Figure 3b shows the change in the concentration values of the second fraction (fraction B) in water. Within the framework of the experiment, the percentage composition of fractions A and B in bottom sediments was calculated (Figure 3c,d). The scenario approach allows us to study the dynamics of changes in the geometry and granulometric composition of the bottom. The software package allows us to model the process of formation of sediments and structures of complex shape, to analyze the transport of suspended matter in a shallow water body, including the estuary area. The developed software package allows us to study the effect of hydrophysical processes on hydrobiocenoses, control the level of water pollution, and also determine the trophic status of a reservoir.

Numerical Implementation of a Hydrophysical Model with Significant Density Gradient in the Aquatic Environment
Consider the results of software implementation of numerical solution of the mathematical model of estuary area hydrodynamics in the form (1)-(4). This model is characterized by of a significant density gradient of the aquatic environment. Define the input data: flow speed is 0.2 m/s; deposition rate is 2.042 mm/s (by Stokes); fresh water density under normal conditions is 1000 kg/m 3 ; suspension density is 2700 kg/m 3 . Let the volume fraction of the suspension is 1/17. Define the calculation area. Consider that its length and width is 50 m; the depth is 2 m. In the calculations, horizontal and vertical steps were used, their values were 0.5 and 0.1 m. The time step was set equal to 0.25 s. The time interval was 5 min.
Analyze the results of a numerical experiment of suspended matter transport for the following scenario: we believe that there is a significant gradient in the density of the aquatic environment in the area under consideration. The results obtained are shown in Figures 4  and 5, where the density is reflected on the right in the cross section of the computational domain by the xOz plane. Consider that this plane is located in the center, while y = 25 m.  The software package developed on the basis of the proposed mathematical models of hydrodynamics and suspension transport has a fairly wide range of applications. The developed software allows us to analyze the transport of suspensions lighter than water, and with appropriate parameterization allows us to predict the spatial change in the concentrations of heavy impurities in a shallow water body.

Discussion
After analyzing the existing models and methods intended for modeling hydrodynamic processes and suspension transfer, it was found that the developed complex of 4D mathematical models allows us to increase the accuracy of modeling the processes under study. Many existing models are developed to simulate hydrodynamic processes in deep water bodies and cannot be used for coastal systems that are characterized by a large difference in depths, salinity, significant influence of winds and river runoff. The proposed model of hydrodynamics is more accurate than the known ones, since it considers the bathymetry of the bottom surface, surge phenomena. When constructing a mathematical model of the hydrophysics of the estuary area, important factors were considered that affect the nature of the researched processes, including wind stresses and friction on the bottom, turbulent exchange, the Coriolis force, evaporation, etc.
For calculations, schemes were used that consider the filling of the cells [28]. The method of volume of fluid (VOF) of rectangular cells with a material medium (liquid) was applied in order to improve the accuracy of calculations.
An estimate of the accuracy of solving the hydrophysics problem of the estuary area for different computational grids (for different discretization parameters) is obtained. As a result, it was found that the relative error of the solution with stepwise approximation of the boundaries can reach up to 70%. The development of the proposed numerical method for solving the problem described in the study allows us to reduce the calculation error to 6%. Thickening the computational grid by a factor of 2-8 horizontally and vertically is less efficient, because does not provide such a result. If the central-difference scheme is used to construct a discrete analogue of the considered problem of impurity transport in a shallow water body, researchers have the problem of loss of accuracy if the grid Péclet number takes large values. The problem can be solved by grid thickening, which is often not efficient enough and significantly increases the computational work. Let us take a simple example. In the numerical implementation of the three-dimensional problem of diffusion-convection, a decrease in the Péclet number by 2 times entails the need to reduce the steps in spatial variables by 2 times, and by 4 times for the time step. It is easy to determine that this will lead to an increase in labor intensity by 32 times. On the other hand, the indicated problem can be solved by developing a new class of difference schemes. In the framework of the study, a modified difference scheme is used to discretize the proposed model of hydrophysics of the mouth area in the case of large values of the grid Péclet number. The scheme takes into account the cell occupancy function [18]. It is based on the Upwind Leapfrog and Standard Leapfrog schemes known in the literature, it is their linear combination of a special type. The novelty of the developed scheme is that it includes weight coefficients calculated by minimizing the approximation error. This difference scheme, when solving diffusion-convection problems, practically does not have grid viscosity and, as a consequence, more accurately describe the behavior of the solution in the case of large grid Péclet numbers, and also preserves the smoothness of the solution at the interface when solving hydrodynamic problems with a complex shape of the boundary surface (no oscillations associated with the stepwise approximation of the boundaries). The use of these schemes in solving hydrodynamic allows us to describe more accurately the dynamics of the aquatic environment in estuary areas, then the classical schemes.
The software package is implemented on the basis of a difference scheme with the optimal value of the weight parameter, which made it possible to increase the time step in comparison with the classical approaches.

Conclusions
The paper proposes: a three-dimensional hydrodynamics model of estuary areas for calculating 3D fields of the water flow velocity vector, as well as a model of the transport of particles of a multicomponent suspension. The approximation of the continuous problem of calculating the water bodies velocity field in terms of spatial variables is based on the balance method taking into account the occupancy factors of the calculated cells, which allows us to consider the complex geometry of the coastline and thereby increase the modeling accuracy. The approximation of the model of the transport of suspended particles based on schemes of splitting into a two-dimensional and one-dimensional problems. The discretization of the developed mathematical models based on the difference scheme obtained as a result of a linear combination of the Upwind and the Standard Leapfrog difference schemes with weight coefficients obtained from the condition of minimizing the order of approximation error. For the numerical implementation of the considered hydrophysics models, which allows us to calculate three-dimensional velocity fields of the water bodies, to study the sedimentation process of suspended particles, including in the case of a multicomponent suspension, and also to predict the dynamics of the water movement process in the estuary area in the case of a significant density gradient of the water medium.