Simulation of Vertical Movements of Seawater in Stratified Reservoirs

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

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 [2]. Temperature stratification affects significantly the distribution of organisms in the water column, the transfer and deposition of impurities harmful to bio-organisms. An increase in the temperature of surface waters causes a violation of vertical water exchange and, accordingly, a decrease in aeration of the deep-water zone, a decrease in solubility and oxygen concentration in water. Stratification by density, temperature and chemical composition limits the convective rise of biogenic elements, carbon dioxide and products of incomplete oxidation of organic substances entering the hypolimnion (cold, salty, dense layers of water) as a result of sedimentation of seston into the surface layers of water. From the beginning of stratification until its end, the surface layer is depleted, and the hypolimnion, on the contrary, is enriched with these substances. As a result, physico-chemical stratification causes an uneven distribution of a number of biologically significant substances in depth and is the reason for the self-organization of a complex structure of ecological niches [3]. Mathematical modeling of mechanical, chemical and biological processes occurring in aquatic ecosystems is urgent, it is associated with problems of ecology and life safety of the population of coastal territories.

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 [4]. According to V.I. Vernadsky, one of the most important manifestations of life is the gas exchange of organisms and the environment, mainly respiration processes based on oxygen consumption [5]. Some outstanding Russian scientists, including S.V. Bruevich, whose works were devoted to the development of analytical research methods, the formulation of the basics of hydro- and bio-hydrochemistry, were engaged in the study of hydrobiological processes of reservoirs. G.G. Matishov and V.G. Ilyichev actively study the conditions of optimal exploitation of water resources, develop models of transport of pollutants in water bodies, and investigate the assessment of their impact on the bioresources of the aquatic environment [6].

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 [7][8], it is shown how to effectively use a linear combination of the Upwind and Standard Leapfrog difference schemes with the optimal values of weight coefficients for the approximation of the transport equation. The efficiency of these methods is achieved by optimization of the approximation error of discrete-continuum model, the exact solution of the mass transfer with constant speed. Studies have shown that this approach also extends to hydrodynamic models with a variable (sign-alternating) velocity without the effect of grid viscosity. Another positive quality of such approximations is that they can be used to simulate complex flow structure, e.g., vortex. Currently, numerous researchers use such a scheme for the simulation of turbulent flows. Members of leading foreign research organizations, such as Stanford University, Imperial College London, etc., as well as members of the Institute of Computational Mathematics of the Russian Academy of Sciences E.M. Volodin, A.V. Glazunov, A.S. Gritsun, N.G. Yakovlev and others [9] published works in which mathematical modeling of climatic changes, hydrodynamic and atmospheric processes and phenomena are carried out on the basis of eddy-resolving schemes. The quasi-hydrodynamic approximation of a continuous medium enables to further reduce the time step requirement and increase the spatial resolution of the model with limited computer memory. In practice, a small term proportional to the second derivative in time or density is added to the system of Navier-Stokes and continuity equations. This approach makes it possible to smooth out non-physical fluctuations in mass density and momentum, as well as total energy, brought along the spatial grid faster than the speed of sound.

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 [1][2]:

,

,

,

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 [2]:

(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 [10] defines density of seawater as a function of salinity , temperature and hydrostatic pressure p, which has the form:

(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 [11][5]:

(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 [3], which causes an active consumption of resources of computing devices and developers who create algorithms and programs that numerically implement them.

The problem of setting boundary conditions and their coordinated assignment at alternating speeds was solved by the method of filling control cells [7]. Taking into account the maximum occupancy and expressing its functional dependence on the node number makes it possible to increase the accuracy of approximation of boundary conditions.

Investigations on the theoretical and experimental selection of approximation methods and schemes led to the choice of splitting methods [13] and regularization. The continuity equations were regularized by B. N. Chetverushkin. The given equation was subjected to a difference approximation according to the FTCS scheme. The balance (transfer) equations of momentum, salinity and total energy were replaced by explicit equations obtained by linear combination of various approximations of the transfer operator (leapfrog and Upwind Leapfrog) [12][15]. According to the splitting method, the transfer of momentum (and velocity), salinity, and total energy was approximated in a fractional step. The water velocity variation, which determined the intensity of transport and mass flows, was introduced at the second fractional step by solving the wave equation. Denote , , , , , , . In the split form [16], after the introduction of a regularizer into the continuity equation and semi-discretization on grid , approximation of the partial derivative in time by the FTCS scheme, the first two equations of system (1) can be written as:

, (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 [13]. Summands (8) determining the wave properties during the propagation of the pulse (the right side of the equation) are approximated by the central differences in space:

, (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) [7].

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 [7]; characterizes the occupancy of the area ; — , — .

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:

• building of the projected pulse changing according to the first equation (3);
• approximate calculation of the function of the spatial distribution of pressure as a function of density and temperature from the previous time layer;
• assessment of the density changes according to the second equation (3);
• calculation of the pressure gradient from the new values of density and temperature, correction of the momentum distribution according to the third equation (3);
• finding a new distribution of the total energy and temperature according to the approximated third equation of the system (1).

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

• binary-numeric masks that predetermine the switching of the difference scheme template with a change in the sign of the velocity;
• variable shifts of the indices of neighboring nodes, which make it possible to write the solution to equation (12) for boundary and internal nodes by a single system of computational operations.

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

The matrix equation with a band tridiagonal matrix A (13) is solved by sweeping [18]. For the initial verification of the software implementation, a simplified model of the state of salt water in the form of P. S. Lineikin’s equation was used.

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:

• execution of one step on the transfer equation according to formulas (10);
• calculation of seawater density (unesco_urs);
• calculation of the pressure of seawater located in the gravity field at a given temperature, salinity and density according to empirical equations of state (rhoTS2P);
• calculation of the speed of sound depending on temperature, salinity and pressure according to the standard and the updated UNESCO formula (speed_of_sound);
• estimation of the fluid viscosity when its volumes move under the influence of pressure forces at all points of the spatial grid (ForceOfFriction);
• cyclic variation of state and time variables during vertical movement of salt water (aqua_process);
• initial setting of constant values characterizing the fluid, initial and boundary conditions of its global hydrophysical equilibrium (start, set_parameters);
• formation of band matrix approximating the equation containing pressure (func5);
• solution to the matrix equation by the sweep method (run_sweep_shuttle);
• temperature « total energy conversion (TFromE, EFromT);
• calculation of total energy balance (TotalPower);
• solution to the diffusion-convection-reaction problem, including the approximation of a combination of Leapfrog and Upwind Leapfrog difference schemes (ADR_solver);
• estimation of gradients and time derivatives with respect to central and directional differences, taking into account the change in the sign of the velocity and the pattern of the difference scheme (diff123).

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.

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.