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

Abstract

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 fields of the aquatic environment movement velocity vector, and equation for calculating the variable density. The approximation of the equations for calculating the velocity field by spatial variables is based on the splitting schemes for physical processes with fluid 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 significantly improve the accuracy of modeling suspended particle transport and consider the factors influencing the studied processes.

1. 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 the flow of the Don River [2]. The development and research of basic spatially heterogeneous pollutant transport mathematical models of a continuous type for different types of pollutant sources, including areal, point and linear, in terms of action—instantaneous and continuous, are described in the article by Berlyand M.E. [3].

Despite a significant number of publications, many effects related to spatial heterogeneity of the environment, interspecies interaction, hydrodynamic characteristics of the environment, temperature and oxygen conditions, salinity, and other characteristics were not included in the models of the hydrobiological processes of the reservoir. Accounting for these effects can make a significant contribution to improving the accuracy of forecasts of dangerous natural and human-made phenomena. The widely used models of hydrodynamics and biological kinetics are often subjected to processes of linearization, simplification, and idealization. At the same time, some of the factors that can affect the spatial distribution of the researched processes are neglected, which, in turn, can affect the accuracy of mathematical modeling of the researched objects and phenomena. Simplification of the parametrization of mathematical models, the use of ordinary differential equations in the formulation of the Cauchy problem, a point description of the researched processes are fundamentally incorrect in systems in which the homogeneity of space is significantly violated, which was noted in the Svirezhev Yu.M. and Logofet D.O. work [4].

The results of an experimental study of the process of solid particles gravitational settling in a wide range of Reynolds numbers are presented in the works [5,6]. The authors of these works compared the obtained experimental data with the data on the particle settling rate calculated for based on the Stokes formula at different concentrations of particles, liquid density and Reynolds number. In this case, the results of experimental studies obtained at low concentrations of particles (less than 7.5%) are consistent with the theoretical data obtained by the Stokes formula. As the concentration of particles increases, contradictions are observed between the experimental and theoretical results, which can be explained by the ability of particles to coagulate; therefore, when clustered, the particles settle at a higher rate. In this case, the greater the concentration of particles, the more intense the process of coagulation, and, consequently, the rate of deposition increases.

In the work of Tishkin V.F., Gasilov V.A., and others, modern methods of mathematical modeling are used to study turbulent mixing and the development of hydrodynamic instabilities [7]. The calculation results obtained using qualitative models of hydrodynamic processes are widely used to study hydrophysical processes, which significantly improves the accuracy of their predictive modeling. Klaven A.B. in their work carried out experimental studies on the basis of hydraulic modeling of the channel process and the movement of river flows [8]. Cea L., Vazquez-Cendon M.E. engaged in mathematical modeling of the convective flow movement when calculating the transport of dissolved substances in it. The shallow water equation was used in this case [9]. The work of Li S., Duffy C.J. is devoted to the mathematical modeling of the hydrodynamics of shallow water and the transport of contaminants [10]. In this case, unstructured grids were used to discretize the continuous problem. Kang S., Sotiropoulos F. developed methods and tools for numerical modeling of three-dimensional turbulent free surface flow [11]. The study of the effect of turbulence on the transfer of individual particles in the form of a load in a river bed with a gravel layer was the goal of the work of scientists Paiement-Paradis G., Marquis G., Roy A. Paiement-Paradis G. [12]. Pu J.H., Shao S.D., Huang Y.F. in their work [13] described the numerical and experimental studies of turbulence in shallow flows with an open channel. Hou J. with co-authors have developed a conservative mathematical model of hydrodynamics, numerically implemented on unstructured grids for shallow water flows, taking into account important factors affecting the nature of the studied processes, including wind stress, evaporation, complex bathymetry of the reservoir [14]. The numerical experiments on simulation of particle-laden flows are described in [15]. Lambert B., Weynans L., Bergmann M. propose the Navier–Stokes equations for incompressible flow and partitioned volume penalization-discrete element method solver for ellipsoidal particle motion. Unstructured Cartesian grid three-dimensional hydrodynamic model coupled with a laterally 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 diffusion-convection 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.

2. Materials and Methods

2.1. 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.

2.1.1. Hydrodynamics Model

The considered mathematical model of the hydrodynamics of the estuary area includes [22]:

• equations of motion (Navier–Stokes equations):

  $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial x}+\frac{\partial }{\partial x}\left(\mu \frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial u}{\partial z}\right),$$

  $$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial y}+\frac{\partial }{\partial x}\left(\mu \frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial v}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial v}{\partial z}\right),$$

  (1)

  $$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial z}+\frac{\partial }{\partial x}\left(\mu \frac{\partial w}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial w}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial w}{\partial z}\right)+g,$$
• continuity equation for variable density:

  $$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0,$$

  (2)

  where $\mathbf{V}=\left\{u,v,w\right\}$ is water flow velocity vector for the studied shallow water body [m/s], $\rho$ is the aquatic environment density [kg/m${}^3$], P is the pressure [Pa], g is the acceleration of gravity [m/s${}^2$], $\mu$ and $\nu$ are components of the turbulent exchange coefficient (horizontal and vertical, respectively) [m${}^2$/s], and t is time [s], $x,y$ and z are the values of distances in spatial coordinate directions [m].

The system (1) and (2) is considered under the following boundary conditions:

• at the entrance: $u=u_0,v=v_0,\left(\partial P/\partial \mathbf{n}\right)=0,\left(\partial \mathbf{V}/\partial \mathbf{n}\right)=0$;
• coastal zone and bottom: $\rho \mu \left(\partial u/\partial \mathbf{n}\right)=-{\tau }_x$, $\rho \mu \left(\partial v/\partial \mathbf{n}\right)=-{\tau }_y$, $\left(\mathbf{V},\mathbf{n}\right)=0$, $\left(\partial P/\partial \mathbf{n}\right)=0$;
• upper bound: $\rho \mu \left(\partial u/\partial \mathbf{n}\right)=-{\tau }_x$, $\rho \mu \left(\partial v/\partial \mathbf{n}\right)=-{\tau }_y$, $w=-\omega -{\left(\rho g\right)}^{-1}\left(\partial P/\partial t\right)$, $\left(\partial P/\partial \mathbf{n}\right)=0$.
  where $\mathbf{n}$ is the normal vector directed inside the computational domain, $\omega$ is the evaporation rate of liquid, and ${\tau }_x$ and ${\tau }_y$ are components of the tangential stress.

On the free surface, we set the tangential stress as: $\left\{{\tau }_x,{\tau }_y\right\}={\rho }_{a}C{d}_s\left|\mathbf{w}\right|\left\{w_x,w_y\right\}$, where $C{d}_s$ = 0.0026, where $\mathbf{w}$ is a vector of the wind speed relative to the water [m/s]; ${\rho }_a$ is the density of the atmosphere [kg/m${}^3$]; $C{d}_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: $\left\{{\tau }_x,{\tau }_y\right\}=\rho C{d}_b\left|\mathbf{V}\right|\left\{u,v\right\}$, $C{d}_b=g{n}^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\in \left[0.025,0.2\right]$.

2.1.2. Model of Transport of Multicomponent Suspended Matters

Transport of suspended multicomponent particles is described using the diffusion-convection equation of the following form:

$$\begin{array}{c}\hfill \frac{\partial c_r}{\partial t}+\frac{\partial \left(u{c}_r\right)}{\partial x}+\frac{\partial \left(v{c}_r\right)}{\partial y}+\frac{\partial \left(\left(w+w_{s,r}\right)c_r\right)}{\partial z}=\\ \hfill =\frac{\partial }{\partial x}\left(\mu \frac{\partial c_r}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial c_r}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial c_r}{\partial z}\right)+F_r,\end{array}$$

(3)

where $c_r$, $w_{s,r}$ are the concentration and the speed of sedimentation of r-th fraction of suspension, respectively, [mg/L], [m/s]; $F_r$ is a function that sets the intensity of the distribution of sources of r-th fraction of suspension [mg/(L·s)].

Equation (3) is considered under the following boundary conditions:

• on the free surface: $\left(\partial c_r/\partial z\right)=0;$
• near the bottom surface: $\nu \left(\partial c_r/\partial z\right)=-w_{s,r}{c}_r$;
• on the lateral boundary: $\left(\partial c_r/\partial \mathbf{n}\right)=0$, if $\left(\mathbf{V},\mathbf{n}\right)\ge 0$ and $\nu \left(\partial c_r/\partial \mathbf{n}\right)=\left(\mathbf{V},\mathbf{n}\right)c_r$, if $\left(\mathbf{V},\mathbf{n}\right)<0$.

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):

$$w_{s,r}=\frac{g{d}_r^2\left({\rho }_{V,r}-\rho \right)}{18\eta },$$

where $d_r$ is the diameter of the deposited r-th fraction particle [m], ${\rho }_{V,r}$ is the density of the deposited r-th fraction particle [kg/m${}^3$], $\rho$ is the density of the medium [kg/m${}^3$], $\eta$ 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 $\eta$, 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\approx 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:

$$w_{s,r}=\sqrt{\frac{4g{d}_r\left({\rho }_{V,r}-\rho \right)}{3\rho \zeta }},$$

where $\zeta$ is the medium hydraulic resistance coefficient. At $2<Re<500$, the value of the resistance coefficient $\zeta =18.5/R{e}^{0.6}$, and in the case of $500<Re<1500$, the resistance coefficient is $\zeta =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.

2.1.3. Equation of Variable Density

The standard formula is used to calculate the density of the aquatic environment:

$$\rho =\left(1-\sum _{r=1}^R{V}_r\right){\rho }_0+\sum _{r=1}^R{V}_r{\rho }_{V,r},c_r=V_r{\rho }_{V,r},$$

(4)

where $V_r$, ${\rho }_{V,r}$ are the volume fraction and density of r-th fraction of suspension, ${\rho }_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\left(V\right)=0$. The solution of the system of Equation (1) is found taking into account Equation (2).

2.2. Approximation of Hydrodynamics Equations

Introduce a uniform grid for approximation of the 3D hydrodynamics mathematical model (1)–(2):

$$\begin{array}{c}\hfill {\overline{w}}_h=\left\{t^n=n\tau ,x_i=i{h}_x,y_j=j{h}_y,z_k=k{h}_z;n=\overline{0,N_t},i=\overline{0,N_x},j=\overline{0,N_y},k=\overline{0,N_z};\right.\\ \hfill \left.N_t\tau =T,N_x{h}_x=l_x,N_y{h}_y=l_y,N_z{h}_z=l_z\right\},\end{array}$$

where $\tau$ 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 diffusion-convection type selected in the process of splitting. For calculations on the intermediate time layer, we will use the following equations:

$$\begin{array}{c}\hfill \frac{\tilde{u}-u}{\tau }+u\frac{\partial \overline{u}}{\partial x}+v\frac{\partial \overline{u}}{\partial y}+w\frac{\partial \overline{u}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{u}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{u}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{u}}{\partial z}\right),\\ \hfill \frac{\tilde{v}-v}{\tau }+u\frac{\partial \overline{v}}{\partial x}+v\frac{\partial \overline{v}}{\partial y}+w\frac{\partial \overline{v}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{v}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{v}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{v}}{\partial z}\right),\\ \hfill \frac{\tilde{w}-w}{\tau }+u\frac{\partial \overline{w}}{\partial x}+v\frac{\partial \overline{w}}{\partial y}+w\frac{\partial \overline{w}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{w}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{w}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{w}}{\partial z}\right)+g,\end{array}$$

(5)

where $\mathbf{V}=\left\{u,v,w\right\}$ is the velocity vector with components in the previous time layer. In Equation (5), we introduced the notation $\tilde{u},\tilde{v},\tilde{w}$ are the components of the velocity vector in the intermediate time layer; $\overline{u}=\sigma \tilde{u}+\left(1-\sigma \right)u$, $\sigma \in \left[0,1\right]$ is the weight of the scheme.

To solve the second subproblem (calculation of pressure), the equation is used:

$$\frac{{\partial }^{2}P}{{\partial x}^2}+\frac{{\partial }^{2}P}{{\partial y}^2}+\frac{{\partial }^{2}P}{{\partial z}^2}=\frac{\widehat{\rho }-\rho }{{\tau }^2}+\frac{1}{\tau }\left(\frac{\partial \left(\widehat{\rho }\tilde{u}\right)}{\partial x}+\frac{\partial \left(\widehat{\rho }\tilde{v}\right)}{\partial y}+\frac{\partial \left(\widehat{\rho }\tilde{w}\right)}{\partial z}\right),$$

(6)

where $\widehat{\rho }$ and $\rho$ 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:

$$\frac{\widehat{u}-\tilde{u}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial P}{\partial x},\frac{\widehat{v}-\tilde{v}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial P}{\partial y},\frac{\widehat{w}-\tilde{w}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial P}{\partial z},$$

(7)

where $\widehat{u}$, $\widehat{v}$, $\widehat{w}$ 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:

$$o_{i,j,k}=\frac{P_{i,j,k}+P_{i-1,j,k}+P_{i,j-1,k}+P_{i-1,j-1,k}}{4\rho g{h}_z}.$$

(8)

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:

$$P(x,y,z,t)=p(x,y,z,t)+{\rho }_{0}gz,$$

(9)

where ${\rho }_{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:

$$\begin{array}{c}\hfill \frac{\tilde{u}-u}{\tau }+u\frac{\partial \overline{u}}{\partial x}+v\frac{\partial \overline{u}}{\partial y}+w\frac{\partial \overline{u}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{u}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{u}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{u}}{\partial z}\right),\\ \hfill \frac{\tilde{v}-v}{\tau }+u\frac{\partial \overline{v}}{\partial x}+v\frac{\partial \overline{v}}{\partial y}+w\frac{\partial \overline{v}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{v}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{v}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{v}}{\partial z}\right),\\ \hfill \frac{\tilde{w}-w}{\tau }+u\frac{\partial \overline{w}}{\partial x}+v\frac{\partial \overline{w}}{\partial y}+w\frac{\partial \overline{w}}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{w}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{w}}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial \overline{w}}{\partial z}\right)+\\ \hfill +g\left(\frac{{\rho }_0}{\rho }-1\right),\end{array}$$

(10)

$$\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2}=\frac{\widehat{\rho }-\rho }{{\tau }^2}+\frac{1}{\tau }\left(\frac{\partial \left(\widehat{\rho }\tilde{u}\right)}{\partial x}+\frac{\partial \left(\widehat{\rho }\tilde{v}\right)}{\partial y}+\frac{\partial \left(\widehat{\rho }\tilde{w}\right)}{\partial z}\right),$$

(11)

and

$$\frac{\widehat{u}-\tilde{u}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial p}{\partial x},\frac{\widehat{v}-\tilde{v}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial p}{\partial y},\frac{\widehat{w}-\tilde{w}}{\tau }=-\frac{1}{\widehat{\rho }}\frac{\partial p}{\partial z},$$

(12)

Note that the term $g\left({\rho }_0/\rho -1\right)$ 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.

2.3. Approximation of the Suspended Particles Transport Model

We approximate by an equation of the diffusion-convection type for the three- dimensional case:

$$\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial x}+v\frac{\partial c}{\partial y}+w\frac{\partial c}{\partial z}=\frac{\partial }{\partial x}\left(\mu \frac{\partial c}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial c}{\partial y}\right)+\frac{\partial }{\partial z}\left(\nu \frac{\partial c}{\partial z}\right).$$

(13)

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) ${\omega }_{\tau }$ was used for calculations. The time step $\tau$ was used for discretization:

$${\omega }_{\tau }=\left\{t_n=n\tau ,n=\overline{0,N_t},N_t\tau =T\right\}.$$

For Equation (13), the splitting schemes into a two-dimensional and one-dimensional problem are used [27,28]:

$$\frac{c^{n+1/2}-c^n}{\tau }+u\frac{\partial \overline{c}}{\partial x}+v\frac{\partial \overline{c}}{\partial y}=\frac{\partial }{\partial x}\left(\mu \frac{\partial \overline{c}}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \frac{\partial \overline{c}}{\partial y}\right),$$

(14)

$$\frac{c^{n+1}-c^{n+1/2}}{\tau }+w\frac{\partial \tilde{c}}{\partial z}=\frac{\partial }{\partial z}\left(\nu \frac{\partial \tilde{c}}{\partial z}\right),$$

(15)

where $c^n$, $c^{n+1/2}$, and $c^{n+1}$ are the values of the concentration fields of a one-component suspension on the current time layer $t^n$, on the intermediate time layer $t^{n+1/2}$, and on the next time layer $t^{n+1}$, respectively, $\overline{c}$ and $\tilde{c}$ are the values of the concentration field of a one-component suspension on a certain intermediate layer located on the segment $t\in \left[t^n,t^{n+1/2}\right]$, $\overline{c}={\sigma }_{xy}{c}^{n+1/2}+\left(1-{\sigma }_{xy}\right)c^n$, $\tilde{c}={\sigma }_z{c}^{n+1}+\left(1-{\sigma }_z\right)c^{n+1/2}$, where ${\sigma }_{xy}$ and ${\sigma }_z$ are the weights of the scheme in problems (14) and (15), respectively.

Solve a model problem of the form (10). When organizing calculations, a grid is used in the form:

$$w_h=\left\{x_i=i{h}_x,y_j=j{h}_y;i=\overline{0,N_x},j=\overline{0,N_y};N_t\tau =T,N_x{h}_x=l_x,N_y{h}_y=l_y\right\},$$

where $h_x$, $h_y$, $\tau$ 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$: $\left\{x\in (x_{i-1/2},x_{i+1/2}),y\in (y_{j-1/2},y_{j+1/2})\right\}$, $q_1-D_1:D_0\cap \left\{x\ge x_i\right\}$, $q_2-D_2:D_0\cap \left\{x\le x_i\right\}$, $q_3-D_3:D_0\cap \left\{y\ge y_j\right\}$, $q_4-D_4:D_0\cap \left\{y\le y_j\right\}$.

Denote ${\mathsf{\Omega }}_m$$m=\overline{0,4}$ the parts of the areas $D_m$ that are filled. Calculate the coefficients $q_m$ based on the works [29]:

$${\left(q_m\right)}_{i,j}=\frac{S_{{\mathsf{\Omega }}_m}}{S_{D_m}},{\left(q_0\right)}_{i,j}=\frac{o_{i,j}+o_{i+1,j}+o_{i+1,j+1}+o_{i,j+1}}{4},{\left(q_1\right)}_{i,j}=\frac{o_{i+1,j}+o_{i+1,j+1}}{2},$$

$${\left(q_2\right)}_{i,j}=\frac{o_{i,j}+o_{i,j+1}}{2},{\left(q_3\right)}_{i,j}=\frac{o_{i+1,j+1}+o_{i,j+1}}{2},{\left(q_4\right)}_{i,j}=\frac{o_{i,j}+o_{i+1,j}}{2}.$$

Let the boundary condition have the form $\frac{\partial c}{\partial \mathbf{n}}={\alpha }_{\mathbf{n}}c+{\beta }_{\mathbf{n}}$. Discretization the convection and diffusion operators is calculated on the basis of [29]:

$$\begin{array}{c}\hfill {\left(q_0\right)}_{i,j}u\frac{\partial c}{\partial x}\simeq {\left(q_1\right)}_{i,j}{u}_{i+1/2,j}\frac{c_{i+1,j}-c_{i,j}}{2h_x}+{\left(q_2\right)}_{i,j}{u}_{i-1/2,j}\frac{c_{i,j}-c_{i-1,j}}{2h_x},\\ \hfill {\left(q_0\right)}_{i,j}\frac{\partial }{\partial x}\left(\mu \frac{\partial c}{\partial x}\right)\simeq {\left(q_1\right)}_{i,j}{\mu }_{i+1/2,j}\frac{c_{i+1,j}-c_{i,j}}{h_x^2}-{\left(q_2\right)}_{i,j}{\mu }_{i-1/2,j}\frac{c_{i,j}-c_{i-1,j}}{h_x^2}-\\ \hfill -\left|{\left(q_1\right)}_{i,j}-{\left(q_2\right)}_{i,j}\right|{\mu }_{i,j}\frac{{\alpha }_x{c}_{i,j}+{\beta }_x}{h_x}.\end{array}$$

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].

2.4. Determination of the Optimal Scheme Weight

Consider the one-dimensional diffusion-convection equation:

$$\frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=\mu \frac{{\partial }^{2}q}{{\partial x}^2},u=const,\mu =const,0<x<L,t>0,$$

(16)

with an initial distribution $q(x,0)=q_0\left(x\right)$ and boundary conditions $q(0,t)=0$, $q(L,t)=0$, $t>0$.

Let us find a solution to the problem of the form (16) in the class of functions$q(x,t)\in C^2\left(\left(0,L\right)\right)\cap C\left(\left[0,L\right]\right)\cap C^1\left(\left(0,\infty \right)\right)\cap C\left(\left[0,\infty \right)\right)$.

For the function $q(x,t)$, set the sum of the Fourier series (finite, in complex form) according to the formula:

$$q(x,t)=\sum _{m=-N}^N{C}_m\left(t\right)\exp \left(j\omega mx\right),$$

(17)

where $\omega =\frac{\pi }{L}$, m is the number of the harmonic, $C_m\left(t\right)=\frac{2}{L}\underset{0}{\overset{L}{\int }}q(x,t)\exp \left(-j\omega mx\right)dx$ is the complex amplitude of m-th harmonic, $j=\sqrt{-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: ${\omega }_h=\left\{x_i=ih;i=\overline{0,N};Nh=L\right\}$. 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:

$$\begin{array}{c}\hfill \frac{\partial q}{\partial t}=-{\left(\mathsf{\Lambda }q\right)}_i,\mathrm{define}{\left(\mathsf{\Lambda }q\right)}_i=u\frac{q_{i+1}+q_{i-1}}{2h}+\mu \frac{-q_{i+1}+2q_i-q_{i-1}}{h^2},\\ \hfill q_i=\sum _{m=-N}^N{C}_m\exp \left(\frac{j\pi mi}{N}\right).\end{array}$$

Note that the eigenvalues of the operator ${\left(\mathsf{\Lambda }q\right)}_i$ can be found by the formula:

$${\lambda }_m=-\left(2\mu \frac{1-\cos \left(\frac{\pi m}{N}\right)}{h^2}+ju\frac{\sin \left(\frac{\pi m}{N}\right)}{h}\right).$$

(18)

Consider a space-time grid $\omega ={\omega }_h\times {\omega }_{\tau }$ by introducing an additional time grid ${\omega }_{\tau }=\left\{t_n=n\tau ,n=\overline{0,N_t},N_t\tau =T\right\}$, where $\tau$, $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:

$$\frac{C_m^{n+1}-C_m^n}{\tau }=-{\lambda }_m\left(\sigma C_m^{n+1}+\left(1-\sigma \right)C_m^n\right)\mathrm{or}{C}_m^n=\left(1-\frac{{\chi }_m}{1+{\chi }_m\sigma }\right)C_m^{n-1},$$

(19)

where $C_m^{n+\sigma }=\sigma C_m^{n+1}+\left(1-\sigma \right)C_m^n$, $\sigma \in \left[0,1\right]$ is the weight of the scheme, ${\chi }_m={\lambda }_m\tau$.

Let us find the value of the weight $\sigma$ 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\left(t^n\right)$ and the approximate value of the function $C_m^n$, can be expressed in terms of the following function:

$${\psi }_m^n=C_m^n-C_m\left(t^n\right),$$

or

$${\psi }_m^n=\left(1-\frac{{\chi }_m}{1+{\chi }_m\sigma }\right){\psi }_m^{n-1}+C_m\left(t^{n-1}\right)\left(1-\frac{{\chi }_m}{1+{\chi }_m\sigma }-\exp \left(-{\chi }_m\right)\right),$$

where $C_m\left(t^n\right)=C_m\left(t^{n-1}\right)\exp \left(-{\chi }_m\right)$.

The estimate of the error in the numerical solution of the problem (16) based on [18]:

$$\left|{\psi }_m^n\right|\le \underset{k=\overline{0,n-1}}{max}\left|C_m\left(t^k\right)\right|\left|\left(\left(1-\exp \left(-{\chi }_m\right)\right)\sigma +\frac{1-\exp \left(-{\chi }_m\right)-{\chi }_m}{{\chi }_m}\right)\right|.$$

(20)

It is obvious that the approximation error with respect to the time variable on the i-th layer does not exceed ${\displaystyle \underset{m}{max}}\left|{\psi }_m^n\right|$.

Let us introduce the notation:

$$z_{1,m}=1-\exp \left(-{\chi }_m\right),z_{2,m}=\frac{1-\exp \left(-{\chi }_m\right)-{\chi }_m}{{\chi }_m}.$$

Then, taking into account (20) and Parseval’s theorem, we obtain:

$$\begin{array}{c}\hfill {∥{\mathsf{\Psi }}^n∥}^2=\sum _{m=-N}^N{\left|{\psi }_m^n\right|}^2\le \sum _{m=-N}^N\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2{\left|z_{1,m}\sigma +z_{2,m}\right|}^2=\\ \hfill =\sum _{m=-N}^N\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2\left[\left({\left(\mathrm{Re}{z}_{1,m}\right)}^2+{\left(\mathrm{Im}{z}_{1,m}\right)}^2\right){\sigma }^2+\right.\\ \hfill \left.+2\left(\mathrm{Re}{z}_{1,m}\mathrm{Re}{z}_{2,m}+\mathrm{Im}{z}_{1,m}\mathrm{Im}{z}_{2,m}\right)\sigma +{\left(\mathrm{Re}{z}_{2,m}\right)}^2+{\left(\mathrm{Im}{z}_{2,m}\right)}^2\right].\end{array}$$

Consider the function:

$$\begin{array}{c}\hfill s_m\left({\chi }_m,\sigma \right)=\sum _{m=1}^N\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2\left[\left({\left(\mathrm{Re}{z}_{1,m}\right)}^2+{\left(\mathrm{Im}{z}_{1,m}\right)}^2\right){\sigma }^2+\right.\\ \hfill \left.+2\left(\mathrm{Re}{z}_{1,m}\mathrm{Re}{z}_{2,m}+\mathrm{Im}{z}_{1,m}\mathrm{Im}{z}_{2,m}\right)\sigma +{\left(\mathrm{Re}{z}_{2,m}\right)}^2+{\left(\mathrm{Im}{z}_{2,m}\right)}^2\right].\end{array}$$

To find $\underset{\sigma }{min}\left|s_m\left({\chi }_m,\sigma \right)\right|$, we define the extremum points of the function $s_m\left({\chi }_m,\sigma \right)$ by the variable $\sigma$:

$$\begin{array}{c}\hfill \frac{\partial s_m\left(\chi ,\sigma \right)}{\partial \sigma }=2\sum _{m=1}^N\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2\left[\left({\left(\mathrm{Re}{z}_{1,m}\right)}^2+{\left(\mathrm{Im}{z}_{1,m}\right)}^2\right)\sigma +\right.\\ \hfill +\left.\left(\mathrm{Re}{z}_{1,m}\mathrm{Re}{z}_{2,m}+\mathrm{Im}{z}_{1,m}\mathrm{Im}{z}_{2,m}\right)\right]=0.\end{array}$$

(21)

From (21) obtain:

$$\sigma =-\frac{{\displaystyle \sum _{m=1}^N}\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2\left(\mathrm{Re}{z}_{1,m}\mathrm{Re}{z}_{2,m}+\mathrm{Im}{z}_{1,m}\mathrm{Im}{z}_{2,m}\right)}{{\displaystyle \sum _{m=1}^N}\underset{k=\overline{0,n-1}}{max}{\left|C_m\left(t^k\right)\right|}^2\left({\left(\mathrm{Re}{z}_{1,m}\right)}^2+{\left(\mathrm{Im}{z}_{1,m}\right)}^2\right)}.$$

(22)

The above reasoning allows us to determine the optimal value of the scheme weight $\sigma$ at which the error in the numerical solution of problem (16) is minimal.

3. 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).

3.1. 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\times 100\times 40$ computational nodes is introduced. The parameters of the computational domain are: $50\times 50\times 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.

3.2. 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.

3.3. 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 Figure 4 and Figure 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. Figure 4 and Figure 5 show the values of the suspended matter concentration averaged over depth on the left. The calculated interval was 1 min and 5 min after the start of the computational experiment. Figure 4 and Figure 5 allow to study the hydrophysical processes of the estuary area, including the transport of suspended matter, the vertical sections on the right show the change in the concentration of suspended matter. Consider that the layers of the aquatic environment are stratified. The density of the aquatic environment changes dynamically over time.

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.

4. 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.

5. 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.