Introduction. In the field of computational mathematics, there are many ways to approximate the model of fluid mechanics. Methods and estimates of approximation quality criteria, such as stability and convergence, are developed, while a combination of approaches to constructing economical difference schemes, such as splitting by physical processes, regularization by B. N. Chetverushkin, a linear combination of the Upwind and Standard Leapfrog difference schemes in aggregate has not been implemented and evaluated before. The authors were faced with the task of approximating each part of the hydrodynamic model split by physical processes with the most adequate scheme and further investigating the correctness of this approach.

Materials and Methods. The mathematical model of hydrophysical processes is closed by the empirical equation of the state of salt water. Significant properties were selected, a mathematical model was built. Difference operators approximated differential operators. An algorithm for layer-by-layer modeling of transients was constructed. The algorithm has been implemented in the form of the program, which mainly contains elementwise (massively-parallel) operations.

Results. Mathematical models of hydrodynamic processes in reservoirs were obtained, taking into account three equations of motion in the presence of a density gradient of the aqueous medium when hydrostatic approximation was abandoned. A new method of calculating the pressure field using B. N. Chetverushkin’s regularizers in the continuity equation was tested. A software module for numerical simulation of hydrophysical processes of water movement with different salinity and density was developed. This is open-source software that provides not only the redefinition of empirical dependences (as algebraic functions), but also the connection of external simulating modules to display dependences algorithmically.

Discussion and Conclusion. The developed model of hydrophysics, taking into account the properties of salt water and the dynamic relationship of the mechanical movement of water with salinity, can be used to study the formation of a nonequilibrium distribution of parameters and identify the most stable parameters of the aquatic environment. The model explains the downward movement of oxygen. That will help in the future to estimate the values of the parameters of the aquatic environment, which are difficult to measure directly. It can be used in the procedure of parametric identification of hard-to-measure parameters of the aquatic environment.

Введение. В области вычислительной математики известно множество способов аппроксимации модели механики жидкости. Учеными выработаны методы и оценки критериев качества аппроксимации, таких как устойчивость и сходимость. Комбинация подходов построения экономичных разностных схем, таких как расщепление по физическим процессам, регуляризация по Б. Н. Четверушкину, линейная комбинация разностной схемы «кабаре» и «крест» в совокупности ранее не реализовывалась и не оценивалась. Перед авторами стояла задача аппроксимировать каждую часть расщеплённой по физическим процессам модели гидродинамики наиболее адекватной схемой и далее исследовать корректность данного подхода.

Материалы и методы. Математическая модель гидрофизических процессов замыкается эмпирическим уравнением состояния соленой воды. Выбираются значимые свойства, строится математическая модель. Разностные операторы аппроксимируют дифференциальные операторы. Строится алгоритм послойного моделирования переходных процессов. Алгоритм реализован в виде программы, которая, в основном, содержит поэлементные (массивно параллельные) операции.

Результаты исследования. Получены математические модели гидродинамических процессов в водоемах, учитывающие три уравнения движения при наличии градиента плотности водной среды при отказе от гидростатического приближения. Апробирован новый способ вычисления поля давления с применением регуляризаторов по Б. Н. Четверушкину в уравнении неразрывности. Разработан программный модуль численного моделирования гидрофизических процессов движения воды с различной солёностью и плотностью. Это открытое программное обеспечение, допускающее не только переопределение эмпирических зависимостей (как алгебраических функций), но и подключение внешних моделирующих модулей для отображения зависимостей алгоритмически.

Обсуждение и заключение. Разработанная модель гидрофизики, учитывающая свойства солёной воды и динамическую связь механического движения воды с солёностью, может применяться для изучения формирования неравновесного распределения параметров и идентификации наиболее стабильных параметров водной среды. Модель объясняет нисходящее движение кислорода, что позволит в будущем оценивать величины параметров водной среды, которые сложно измерить непосредственно. Она может быть использована в процедуре параметрической идентификации трудноизмеряемых параметров водной среды.

Introduction. One of the critical tasks related to the ecology and life safety of people living in coastal areas is forecasting and modeling the movement of water in the seas and large regional reservoirs. Moreover, it is important to study the process of transport of substances dissolved in an aqueous medium, taking into account stratification and the dependence of water density on many variable factors. Such predictive modeling is likely to make it possible not only to assess water quality, but also to manage it under conditions of climate change and industrial impacts. A more general and major goal is to link water quality to the number and species diversity of natural hydrobiocenoses living in the hydrosphere. When solving these problems, it is necessary to take into account the hydrodynamic characteristics of the aquatic environment, features of external influencing factors, such as the heterogeneity of temperature distribution, salinity, oxygen saturation of water, the amount of gases dissolved in water, acidity. These parameters are some key parameters of the biological activity of the aquatic ecosystem [

The water area of the reservoir can also be considered as a transport system that transfers oxygen from the atmosphere to bottom sediments. However, there are cases of the formation of chemical gradients of large magnitude at relatively small depth differences — in boundary layers with a relatively small value outside these zones. The reason for the appearance of such sites with the general stratification of water by density in universal gravity on the one hand, and the radiation effect of the sun on the water, causing its heating, on the other. These processes can cause a decrease in the rate of production, destruction and recycling of biogenic elements and bio-organisms until these processes stop, as well as predetermine the biodiversity of hydrobionts in general and species composition in particular [

We should name the outstanding scientists who made a significant contribution to the study of hydrology and oceanology. V.P. Dymnikov was engaged in the study of climate and oceanology, modeling of the atmosphere and ocean. A.S. Monin and M.Y. Belevich investigated the processes of kinematics of the aquatic environment, turbulence and microstructure of the ocean. Soviet scientist A.S. Sudolsky studied the dynamics of waters and coastal processes in various reservoirs in relation to solving the problems of designing, building and operating specific hydraulic structures, and rational economic use of reservoirs in general [

To study the influence of these processes of the aquatic environment, a complex of interrelated mathematical models is being developed. It is based on the use of accurate predictive models and software implementation of economical numerical methods, which provide detailed investigation of the kinematics of the process, cause-effect relationships, and the state of the modeling object. The existing methods and means of predictive modeling of the state of the aquatic environment, taking into account a number of biotic and abiotic factors, including the processes of distribution of oxygen, carbon dioxide, salts, are based on general scientific approaches, simplified mathematical models with low adaptability, lack of modeling of nonlinear dynamic processes characteristic of most aquatic ecosystems, incorrect setting of the boundaries of the computational domain. In some cases, the studies are accompanied by a formal definition of boundary conditions, and give rather rough and approximate simulation results.

When modeling the substance transfer process based on the advection-diffusion equations, a good approximation of advective terms, which are gradients of pressure, mass density and total energy, momentum of motion, is required. The use of standard difference schemes with large estimates of similarity parameters causes loss of calculation accuracy due to an increase in the approximation error and increased restrictions on the time step due to the stability condition of the difference scheme. In the works of A.I. Suhinov, A.E. Chistyakov and others [

Existing universal application software packages (e.g., “Mars3d” software package, Ecointegrator, CHARISMA software package, SALMO complex, CHTDM, CARDINAL software package, packages for modeling various processes of aerohydrodynamics, PHOENICS, FLUENT, GAS DYNAMICS TOOL software packages) do not take into account some properties of the simulated complex systems, thus reducing the accuracy and efficiency of modeling. These properties include spatial heterogeneity of the aquatic environment motion, vortex structures of currents. Mathematical models and algorithms for their numerical implementation do not take into account the probability of a significant change in the depth and density of the aquatic environment, which can cause instability of the numerical solutions obtained. For this reason, such specialized software packages can be used to simulate a limited variety of hydrophysical processes of water systems. Most of the well-known specialized software (ADAM, CAL3QHC, Chensi, TASCflow, ISC–3, PANACHE, REMSAD, UAM-IV, ECOLOGIST, PRISM, VITECON), designed to simulate the process of spreading pollutants, the interaction of hydrobionts, is mainly focused on uniprocessors, represented mainly by personal computers. In such systems, only single component modules of these systems (e.g., ECOSIM and MAQSIP) are scaled to parallel systems. In practice, a small term proportional to the second derivative in time or density is added to the system of Navier-Stokes equations and continuity. This approach makes it possible to smooth out non-physical fluctuations in mass density, momentum and total energy brought along the spatial grid faster than the speed of sound.

Materials and Methods. The success of the development of a mathematical model of hydrophysical and hydrobiological processes depends on the availability and elaboration of test cases and tasks for investigating steadily observed phenomena in the seas, such as vertical mixing and redistribution of salinity and oxygen, halocline and thermocline. To study these phenomena, the paper uses a model of hydrodynamics that takes into account the balance of mass forces and cross-border flows [

,

,

,

where — flow density; — kinetic-energy density rate; — self-energy density rate; — stress tensor ; — mass force; — heat flux density; — specific heat inflow due to radiation. Due to the fact that these phenomena are usually described by a vertical distribution of parameters over depth, it is reasonable to obtain a simplified model that provides rapid identification of unobservable parameters. We select a cylindrical region V with bases on the bottom and surface of the water with a cross-sectional area S. We project the velocities, flows and forces to the vertical direction, assuming that the partial derivatives in x, y of the parameters to the horizontal direction are equal to zero outside the cylinder, i.e., assume horizontal uniformity of the parameters of the aqueous medium. Let us write system (1) in a conservative form so that it enables to determine the mass density (), mechanical momentum () and total energy density (), . In the cylinder, we allocate infinitesimal volume and assume that each such volume, which is V, is affected by the reaction force of the support (bottom of the reservoir), equal to hydrostatic pressure, and a force similar to the friction force caused by the viscosity of the fluid and momentum transfer, not equal to zero with vertical movements of the fluid.

If we neglect horizontal fluid movements and assume that only vertical movements are essential, and the fact that the density of the aqueous medium significantly depends on salinity, accept the simplifications and conventions outlined earlier, then the equations of hydromechanics [1–3] in compact form are presented by a system of partial differential equations [

(1)

where , — the empirical equation of the state of seawater and the equation, closing the system by internal energy, respectively; p — total hydrodynamic pressure; — local density of the aquatic environment; — projection of the velocity vector function on the vertical (z axis is directed upwards from the bottom to the surface); — volume density of internal energy; — pressure of the gas enclosed in an elementary volume between adjacent layers; — volume density of the generalized force (sum of forces) applied to elementary volumes of fluid, in addition to pressure; — salt concentration; — cross-sectional area of the cylindrical selected area in which the most intensive upwelling and salt transport process is assumed; — absolute water temperature; — temperature of the water external to the allocated volume; — thermal conductivity of water; — acceleration of gravity; — coefficient characterizing the intensity of momentum transfer due to viscosity.

The boundary conditions characterizing the property of impermeability of fluid through the rock that makes up the bottom of a water body for system (1) can be written in the form of equalities:

. (2)

Where n — the normal vector directed inside the computational domain.

The International Thermodynamic Equation of Seawater [

(3)

where — average modulus of elasticity; digit 0 corresponds to one standard atmosphere (101,325 Pa).

Under numerical modeling based on (1) – (3), also in demand as (3), we establish an algebraic relationship of pressure with density, temperature and salinity. The quadratic equation has two roots:

(4)

(5)

where , , , — variable parameters whose relation to temperature, salinity and density is determined by the standard model of seawater; , where — water Celsius temperature. The speed of sound in salt water can be expressed by the formula: .

Suppose that the momentum, salinity and heat of a volume of water, bounded by a cylindrical surface, change only when the parameters at the boundary of the domain (Dirichlet boundary conditions) change, in particular, the salinity changes on the surface of the water. The behavior of mechanical parameters — density, momentum — is set by Neumann conditions and is inaccessible to external influence. Regardless of these assumptions, within the boundary of the aqueous medium, the system of equations (1) can be written in the vector form of a transfer-reaction model with a source additive [

(6)

where — vector characterizing the interaction of the flow with the surrounding fluid and the planet; — vector of conservative state variables; — vector of flows acting as feedback and closing the balance equation. From the components of vector common multiplier can be derived .

In practice, to solve the transfer equation of form (6), a difference scheme has proven itself well, whose layered transition operator is obtained by a linear combination of similar transition operators of the Upwind and Standard Leapfrog schemes [12–15]. Taking into account the relationship between mass flows and density changes turned out to be more effective for a quasi-hydrodynamic system regularized according to B.N. Chetverushkin, approximated according to the FTCS scheme (forward-time central-space).

The discretization of the continuum model will be carried out by the integro-interpolation method on a uniform grid , , where — vertical grid step, i — index of the node (control/final volume in terms of the Godunov method) when numbering in space, — grid step in time, — the number of the time layer. We describe a finite-difference approximation of the model. The difference schemes used to approximate the first-order balance equations give relatively acceptable results only with a very small grid step [

The problem of setting boundary conditions and their coordinated assignment at alternating speeds was solved by the method of filling control cells [

Investigations on the theoretical and experimental selection of approximation methods and schemes led to the choice of splitting methods [

, (7)

,

, (8)

(9)

Equations (7) and (9) are a discrete analogue of the momentum transfer and change model — the second equation of system (1). Equation (8) predicts a change in the pressure field taking into account the continuity of the flow and information about internodal (intercellular) and boundary flows. Assuming time constancy at each step of the speed of sound, equation (8) can be approximated by an implicit time difference scheme and solved by sweeping. It was this direction in numerical modeling that was selected for the study.

Equation (7) and the fourth equation (1) are approximated on an equidistant grid by an explicit difference scheme [

, (10)

where — occupancy of the cell to the right of the node with index ; — occupancy of the cell to the left of the node with index ; — degree of occupancy of the control area , located in the neighborhood of the node with index i (, , — index-dependent variables) [

The remaining two summands (8) determine the relationship of the rate of change in density with its flow. (8) includes the ratio of the rate of change in density to the time interval, expressed by the difference with a forward shift:

,

the mass flow is transferred to the right side (8) and approximated by the central differences:

(11)

where , , — degree of occupancy of the areas located in the vicinity of the cells with number [

Balance equations of the mechanical impulse and pressure (9) are approximated by FTCS:

(12)

Let us present the problem of solving the difference approximation of equation (8) in the form of a matrix sweep problem with a vector variable in time on the right side :

.

The approximation of equation (8), defined on a three-point difference template, in the form of a linear system of equations , solvable with respect to pressure (p), has the form:

(13)

A similar band matrix is obtained during a similar approximation of the other two terms (8).

The model can take into account dynamic changes in the flow at the boundaries and other parameters, such as temperature, salinity, oxygen content. Therefore, the simulation of water movement must be performed in a time cycle in layers with the retention of operational information about the parameters on at least two time-adjacent layers of the solution to the grid equation. The algorithm for calculating hydrodynamic parameters on a two-index grid in space and time includes:

When algorithmizing the method for solving a system of split equations, the following are introduced:

Such variable index shifts are used in calculating gradient approximations both in the volume of the continuum model and at the boundaries with the second order of accuracy according to discrete analogues of equations (7) – (9) approximated by FTCS and by a linear combination of Leapfrog and Upwind Leapfrog schemes [12–14]:

(14)

where is characteristic q for transfer of salt () and momentum .

Parameter m is a “switch” of the flow when the sign of the velocity changes in the difference approximations of the pressure gradients and the mass flow:

,

When using such a switch, the approximation of the gradient in solving the transfer problem for momentum and salinity can be written by the formula:

. (15)

The use of the cell occupancy method makes it possible to correctly approximate the boundary conditions when placing water state parameters in computer memory with an array of numerical values, and variable index shifts enable to reduce the volume of symbolic recording of the subroutine for performing layer-by-layer iterative changes of state variables assigned to the nodes of the difference scheme [

The matrix equation with a band tridiagonal matrix A (13) is solved by sweeping [

Research Results. The simulation was performed on the basis of a software package written in MATLAB. Debugging was performed using the GNU Octave interpreter. Operation also presupposes the availability of libraries of this system. The software package consists of 16 functional modules. The selection of this interpreter and the corresponding language was due to the ability to write and test a program that operates with an array of state variables, which is a projection of the desired functions , , on a spatial grid. At this stage of the study, the authors abstracted from the specifics of performing element-by-element operations on arrays of real numbers.

The software system consists of the interconnected subprograms:

The constructed model was used in a test run to estimate the density change under an increase in the salinity of the surface of the aquatic environment by 1% at the 30th second from the start of the simulation and a maximum depth of 1 km. Density transients in a set of cross-sections (function ), spaced from each other, caused by an instantaneous change in salt concentration, are finite in time (Fig. 1), and the density pulse front is greatly weakened when moving in space.

Fig. 1. Graphs of dependence of density on time with a sharp change in salinity on the water surface

Discussion and Conclusion. The presented spatially distributed model does not yet allow predicting a steady redistribution of density and a change in the salinity gradient in numerical calculations, because it does not take into account the buoyancy of less salty water and its change during salinization of the upper layers. Mathematics and software, which greatly simplify the modeling of processes leading to the observed effects of halocline, chemocline and pycnocline, has been tested and debugged on a variety of test tasks (momentum and salt transport, pressure wave propagation). For high-precision modeling of the observed physico-chemical phenomena in the seas, it is required to solve additional tasks of identifying model parameters, taking the data of observations and remote sensing of the Earth as initial information.

The authors declare that there are no conflicts of interest present.