Next Article in Journal
Twisted Edwards Elliptic Curves for Zero-Knowledge Circuits
Next Article in Special Issue
Two-Dimensional-One-Dimensional Alternating Direction Schemes for Coastal Systems Convection-Diffusion Problems
Previous Article in Journal
Implicit Solitary Waves for One of the Generalized Nonlinear Schrödinger Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Modeling of the Phytoplankton Populations Geographic Dynamics for Possible Scenarios of Changes in the Azov Sea Hydrological Regime

by
Alexander Sukhinov
1,
Yulia Belova
1,*,
Alexander Chistyakov
1,
Alexey Beskopylny
2,* and
Besarion Meskhi
3
1
Department of Mathematics and Informatics, Faculty of IT-Systems and Technology, Don State Technical University, 344000 Rostov-on-Don, Russia
2
Department of Transport Systems, Faculty of Roads and Transport Systems, Don State Technical University, 344000 Rostov-on-Don, Russia
3
Department of Life Safety and Environmental Protection, Faculty of Life Safety and Environmental Engineering, Don State Technical University, 344000 Rostov-on-Don, Russia
*
Authors to whom correspondence should be addressed.
Mathematics 2021, 9(23), 3025; https://doi.org/10.3390/math9233025
Submission received: 23 October 2021 / Revised: 21 November 2021 / Accepted: 23 November 2021 / Published: 25 November 2021

Abstract

:
Increased influence of abiotic and anthropogenic factors on the ecological state of coastal systems leads to uncontrollable changes in the overall ecosystem. This paper considers the crucial problem of studying the effect of an increase in the water’s salinity in the Azov Sea and the Taganrog Bay on hydrobiological processes. The main aim of the research is the diagnostic and predictive modeling of the geographic dynamics of the general phytoplankton populations. A mathematical model that describes the dynamics of three types of phytoplankton is proposed, considering the influence of salinity and nutrients on algae development. Discretization is carried out based on a linear combination of Upwind Leapfrog difference schemes and a central difference scheme, which makes it possible to increase the accuracy of solving the biological kinetics problem at large values of the grid Péclet number (Peh > 2). A software package has been developed that implements interrelated models of hydrodynamics and biogeochemical cycles. A modified alternating-triangular method was used to solve large-dimensional systems of linear algebraic equations (SLAE). Based on the scenario approach, several numerical experiments were carried out to simulate the dynamics of the main species of phytoplankton populations at different levels of water salinity in coastal systems. It is shown that with an increase in the salinity of waters, the habitats of phytoplankton populations shift, and marine species invasively replace freshwater species of algae.

1. Introduction

From 2006 to 2020, the average salinity of the Azov Sea changed from 9.4 to 14‰, which could have a beneficial effect on the efficiency of marine fish spawning. However, at the same time, it became the reason for a significant restriction in the distribution of freshwater species of commercial fish. The increase in the salinity of the Azov Sea and the Taganrog Bay becomes more intense due to the lack of water on floodplain spawning grounds for such commercial fish species such as bream and roach. As a result, stocks of bream and roach have been at a consistently low level in recent years. According to the Azov–Black Sea branch of the Russian Federal Research Institute of Fisheries and Oceanography (“AzNIIRKH”), since 2015, the amount of roach stocks has decreased from 8.4 thousand tons to 2.2–2.9 thousand tons and currently does not exceed two thousand tons. In recent years, stocks of bream have averaged 500–600 tons. Among the main reasons for the salinization of the sea are a decrease in the supply of freshwater from the Don and the Kuban Rivers, the main freshwater arteries of the Azov basin, as well as other reasons (e.g., global warming).
The Azov Sea is an example of a coastal system with a high reaction rate to fluctuations in river flow and changes in climatic conditions, which entail a sizeable spatio-temporal variability of hydrophysical and biological parameters of processes. The dynamics of phytoplankton populations are an essential indicator of the reservoir ecosystem state since phytoplankton makes up about 90% of the biomass of the sea and, along with zooplankton, is located at the base of the food pyramid. A long-term study of the Azov Sea phytoplankton revealed roughly 103 species and intraspecific taxa of cyanoprokaryotes (blue-green algae), which account for up to 90% of the total phytoplankton biomass [1]. Phytoplankton biomass in the Azov Sea varies from 3–10 to 70–80 g/m3. The most common types are: Aphanizomenon flos-aquae (L.) Ralfs, Microcystis aeruginosa Kutz. emed Elenk., Anabaena flos-aquae Elenk, and Lyngbya limnetica Lemm. In the last decade, in the plankton of the eastern regions of the Azov Sea and the Taganrog Bay, several species are often recordedt: Planktothrix agardhii (Gom.) Anagn. et Kom, and Nodularia spumigena Mertens [2]. The appearance of these species may be associated with a periodic increase in salinity in the delta region. The mass development of blue-green algae in Taganrog Bay usually begins in July, at a water temperature of 24–29 °C, and reaches a maximum in early autumn. The value of salinity, along with biogenic, temperature, and oxygen regimes, refers to the limiting factors of the development of microalgae in estuarine water bodies. Most cyanoprokaryotes should be classified as freshwater organisms, but the euryhaline content allows them to vegetate intensively far beyond the freshened zone. The mass development of cyanoprokaryotes in the Azov Sea is limited by the 7–10.5‰ isohaline [3]. Freshwater Aphanizomenon flos-aquae can exist in a salinity range from 1.5 to 12‰. However, the growth of its biomass and maximum abundance are noted in the salinity range of 1.5–6.5‰.
From 2014–2015 in the Azov-Black Sea basin, an abnormal frequency of stormy southeastern winds and the advection of salty Black Sea waters into the Azov Sea was observed. During intense storm cyclones in the northeast of the Black Sea, water with a salinity of 15–17‰ was injected into the basin of the Azov Sea. As a result, Azov water with a salinity of about 11–14‰ was displaced towards the Don River runoff into the Taganrog Bay. In summer 2015, there was a decrease in the flow of the Don River into the water area of the bay, which also contributed to an increase in salinity (up to 1.5–2.5‰). With an increase in the concentration of salt in the environment, the cells of blue-green algae, similar to other organisms, experience osmotic and ionic stress and secrete odorants, which cause an unpleasant odor in the water. The cause of the unpleasant smell of water is not the “bloom” itself, but the mass dying off of microalgae in an environment of water salinity abnormal for the seaside (2–5‰).
Many Russian and foreign scientists are studying the effect of changes in the hydrological regime on aquatic organisms in the aquatic environment. Therefore, in [4] the expert system “Lakes of Karelia”, which takes into account the influence of external factors, including salinity and temperature, on the productivity of aquatic ecosystems, is described. Research [5] is devoted to the seasonal, spatial distribution of phytoplankton chlorophyll “a” in the continental shelf ecosystem. According to the review, it was found that changes in the salinity of an aquatic ecosystem can lead to a change in its species composition and the displacement of freshwater species of aquatic organisms by marine ones with an increase in its concentration [6,7]. The salinity and temperature regime is an essential link in the productivity and reproduction of valuable and commercial fish of aquatic biogeocenosis. Therefore, its study based on the methods and means of mathematical modeling is a problem of great national economic importance [8].
Many foreign scientists study the processes of hydrodynamics, the influence of the salinity regime on the biota in coastal systems, whose features are a large difference in depths, a significant effect of wind influences, currents, and river flows. French scientists, together with Russian researchers, have developed a three-dimensional hydrodynamic model (MARS) [9], which describes the processes occurring in marine systems. Additional modules have been developed to simulate sediment dynamics, microbiology, and the contaminants spread, biogeochemical cycles (primary production, development of toxic plankton species, oxygen regime, etc.). Work [10] is devoted to modeling salinity and bottom currents in the shallow Mediterranean lagoon, in which the results of a computational experiment are presented. The influence of water exchange between water bodies with different salinity levels has been investigated. In [11], an ecological model, ECO-MARS3D, of the pelagic ecosystem of the Biscay Bay and the continental English Channel shelf (northeastern Atlantic) was proposed for modeling eutrophication, the presence of contaminants, the level of dissolved oxygen in coastal waters. The model includes inorganic nutrients (NH4, NO3, PO4, Si(OH)4), three main types of phytoplankton (diatoms, nanoflagellates, and dinoflagellates), two main types of zooplankton, particulate matter, and dissolved oxygen.
In work [12], the estuary—located in the joint zone between the Jiaojiang River and ocean—was researched. A two-dimensional mathematical model MIKE 21 [13] was used to ascertain the dynamics of salinity and nutrients in the studied estuary and understand the influence of draining on hydrobionts in this area. The researchers divided the sea area into three parts: nearshore area with low-salt and high eutrophication, medium mixed salinity under nitrogen limitation area, and high-salt under phosphorus limitation area. Modeling allows determining that the change in salinity and concentration of nutrients due to the discharge of a large volume of water during the operation of the sluice of the estuary water area affects the development of aquatic organisms but would not be lethal.
In work [14], the hydrodynamic characteristics and tidal phenomena in estuaries are investigated. The study uses a simplified estuary model to improve the understanding of oscillation and tidal wave reflection behavior in estuaries. An analytical energy approach and a hydrodynamic numerical (HN) model were used to determine the reflection coefficients, and harmonic analysis was used to analyze the simulation data.
Biogeochemical processes in marine ecosystems are modeled by Norwegian scientists together with Russian researchers. Therefore, in work [15], the complex of C-N-P-Si-O-S-Mn-Fe transformation model BROM and the 2–dimensional benthic-pelagic transport model (2DBP) were applied to analyze the impact of salmon farms on the water column and benthic biogeochemistry. Hydrophysical model data for the Hardangerfjord in western Norway were applied for the 2DBP model. The model predicted significant impacts on seafloor biogeochemistry up to 1 km from the fish farm. In addition, studies of the hydrochemical regime of waters along the Spitsbergen Archipelago were carried out during joint Norwegian–Russian expeditions in 2014–2015 [16].
Based on the previous analysis, it can be concluded that studying the influence of changes in the hydrological regime on the dynamics of phytoplankton population development in coastal systems is an urgent task. Therefore, diagnostic and predictive modeling of the main phytoplankton population’s geographic dynamics is the primary goal of this research.

2. Materials and Methods

A three-dimensional mathematical model of biological kinetics has been developed to study the effect of salinity on the species diversity of phytoplankton in the coastal systems of the Azov Sea and Taganrog Bay.

2.1. Mathematical Model of the Phytoplankton Populations Dynamics

A multi-species mathematical model describing the dynamics of phytoplankton populations under changing salinity regimes, taking into account the transformation of nutrient forms and their consumption by algae, is based on a system of non-stationary equations of a convection–diffusion reaction of the parabolic type with nonlinear functions of sources and lower derivatives [17,18]:
c i t + u c i x + v c i y + w c i z = d i v ( k grad c i ) + R c i ,
where c i is the concentration of i-th component, [mg/L]; iM, M = {F1, F2, F3, PO4, POP, DOP, NO3, NO2, NH4, Si}; V = { u , v , w } are the components of water flow velocity vector, [m/s]; k is the turbulent exchange coefficient, [m2/s]; R c i are the chemical-biological sources, [mg/(L∙s)]. The shallow depth is a distinctive feature of the investigated reservoir—the Azov Sea. On average, the depth is 7.4 m, the maximum is 13.5 m. As a result, the sea area in the warm season warms up almost evenly. Therefore, bathymetry is not considered in the model of the dynamics of phytoplankton populations.
In Equation (1) the index i indicates the type of substance (Table 1).
The following equations describe chemical-biological sources:
R F i = C F i ( 1 K F i R ) c F i K F i D c F i K F i E c F i ,   i = 1 , 3 ¯ ,
R P O P = i = 1 3 s P K F i D c F i K P D c P O P K P N c P O P ,
R D O P = i = 1 3 s P K F i E c F i + K P D c P O P K D N c D O P , R P O 4 = i = 1 3 s P C F i ( K F i R 1 ) c F i + K P N c P O P + K D N c D O P ,
R N H 4 = i = 1 3 s N C F i ( K F i R 1 ) f N ( 2 ) ( c N H 4 ) f N ( c N O 3 , c N O 2 , c N H 4 ) c F i + i = 1 3 s N ( K F i D + K F i E ) c F i K 42 c N H 4 ,
R N O 2 = i = 1 3 s N C F i ( K F i R 1 ) f N ( 1 ) ( c N O 3 , c N O 2 , c N H 4 ) f N ( c N O 3 , c N O 2 , c N H 4 ) c N O 2 c N O 2 + c N O 3 c F i + K 42 c N H 4 K 23 c N O 2 ,
R N O 3 = i = 1 3 s N C F i ( K F i R 1 ) f N ( 1 ) ( c N O 3 , c N O 2 , c N H 4 ) f N ( c N O 3 , c N O 2 , c N H 4 ) c N O 3 c N O 2 + c N O 3 c F i + K 23 c N O 2 ,
R S i = s S i C F 3 ( K F 3 R 1 ) c F 3 + s S i K F 3 D c F 3 ,
where K F i R is the specific respiration rate of phytoplankton; K F i D is the specific rate of phytoplankton dying; K F i E is the specific rate of phytoplankton excretion; K P D is the specific speed of autolysis POP; K P N is the phosphatification coefficient POP; K D N is the phosphatification coefficient DOP; K 42 is the specific rate of oxidation of ammonium to nitrites in the process of nitrification; K 23 is the specific rate of oxidation of nitrites to nitrates in the process of nitrification; s P , s N , s S i are the normalization coefficients between the content of N, P, Si in organic matter.
The following expressions determine the growth rate of phytoplankton:
C F 1 , 2 = K N F 1 , 2 f T ( T ) f S ( S ) m i n { f P ( c P O 4 ) , f N ( c N O 3 , c N O 2 , c N H 4 ) } ,
C F 3 = K N F 3 f T ( T ) f S ( S ) m i n { f P ( c P O 4 ) , f N ( c N O 3 , c N O 2 , c N H 4 ) , f S i ( c S i ) }
where K N F is the maximum specific growth rate of phytoplankton.
Functions of the dependence of the growth rate of aquatic organisms on temperature and salinity:
f T ( T ) = exp ( a l { ( T T o p t ) / T o p t } 2 ) ,   l = 1 , 3 ¯ ;   f S ( S ) = exp ( b l { ( S S o p t ) / S o p t } 2 ) ,   l = 2 , 3 ;
f S ( S ) = { k s ,       for       S S o p t , exp ( b 1 { ( S S o p t ) / S o p t } 2 ) ,       for       S > S o p t ,
where k s = 1 , T o p t , S o p t are the temperature and salinity optimal for a given type of aquatic organisms; and a l > 0 , b l > 0 are the coefficients of the width of the range of aquatic organisms tolerance to temperature and salinity, respectively.
The study was carried out under conditions typical for the warm season when the Azov Sea’s water temperature fluctuations are 4 °C. However, to make forecasts for the occurrence of abnormal temperature phenomena, it is necessary to take into account this abiotic factor. For this, a temperature field is connected to model (1), constructed similarly to the salinity field (see Section 3) from satellite data on the water surface temperature. T o p t and S o p t are the same for the entire study area since the values of salinity and temperature are optimal for each phytoplankton species. However, salinity has a more significant influence on the phytoplankton population’s development, which varies from 0‰ at the mouth of the Don River to 12–14‰ at the outlet to the Black Sea.
Functions describing biogen content:
-
for phosphorus f P ( c P O 4 ) = c P O 4 c P O 4 + K P O 4 ,
where K P O 4 is the half-saturation constant of phosphates.
-
for silicon f S i ( c S i ) = c S i c S i + K S i ,
where K S i is the half-saturation constant of silicon.
-
for nitrogen f N ( c N O 3 , c N O 2 , c N H 4 ) = f N ( 1 ) ( c N O 3 , c N O 2 , c N H 4 ) + f N ( 2 ) ( c N H 4 ) ,
f N ( 1 ) ( c N O 3 , c N O 2 , c N H 4 ) = ( c N O 3 + c N O 2 ) exp ( K p s i c N H 4 ) K N O 3 + ( c N O 3 + c N O 2 ) ,   f N ( 2 ) ( c N H 4 ) = c N H 4 K N H 4 + c N H 4
where K N O 3 is the half-saturation constant of nitrates, K N H 4 is the half-saturation constant of ammonium, K p s i is the coefficient of ammonium inhibition.
An initial-boundary value problem is posed in a cylindrical domain G for the system (1). It is assumed that the boundary ∑ of the cylindrical region G is a piecewise-smooth surface and ∑ = ∑H ∪ ∑oσ, where ∑H is the surface of the reservoir bottom ∑o is the undisturbed surface of the water medium, σ is the lateral (cylindrical) surface. Let un is normal with respect to the ∑ component of the water flow velocity vector, n is the vector of the external normal to ∑. The boundary conditions are determined for the concentrations ci:
c i = 0   on   σ , if   u n < 0 ;   c i n = 0   on   σ ,   if   u n 0 ; c i z = 0 ,   on   b ;   c i z = ε i c i   on   H ,
where ε i are the non-negative constants, iM; ε i take into account the sinking of algae to the bottom and their flooding for i ∈ {F1, F2, F3} and take into account the absorption of nutrients by bottom sediments for i ∈ {PO4, POP, DOP, NO3, NO2, NH4, Si}.
It is also necessary to add the initial conditions:
c i ( x , y , z , 0 ) = c i 0 ( x , y , z ) ,   ( x , y , z ) G ¯ ,   t = 0 ,   i M , V ( x , y , z , 0 ) = V 0 ( x , y , z ) ,   T ( x , y , z , 0 ) = T 0 ( x , y , z ) ,   S ( x , y , z , 0 ) = S 0 ( x , y , z ) ,
where c i 0 ( x , y , z ) are the preset functions, fields of water flow vector velocity, salinity, and temperature are input data for this model.

2.2. Numerical Solution of the Diffusion-Convection Problem

Consider the discretization of the problem (1)–(3) in the two-dimensional case using the example of the convection–diffusion reaction equation. In the three-dimensional case, discretization is performed similarly.
c t + u c x + v c y = ( μ c x ) x + ( μ c y ) y + f
with boundary conditions:
c n ( x , y , t ) = α n c + β n ,
where u, v are the velocity vector components, μ is the turbulent exchange coefficient, f is the function describing the intensity and distribution of sources, and αn, βn are the given coefficients.
A uniform grid is introduced:
w h = { t n = n τ ,   x i = i h x ,   y j = j h y ;   n = 0 , N t ¯ ,   i = 0 , N x ¯ ,   j = 0 , N y ¯ ;   N t τ = T ,   N x h x = l x ,   N y h y = l y } ,
where τ is the time step, h x , h y are the space steps, N t is the upper time limit, N x , N y are the space boundaries, l x , l y are the characteristic dimensions of the computational grid.
The splitting schemes in space are used to discretize the homogeneous Equation (4):
c n + 1 / 2 c n τ + u ( c n ) x = ( μ ( c n ) x ) x ,
c n + 1 c n + 1 / 2 τ + v ( c n ) y = ( μ ( c n + 1 / 2 ) y ) y .
When solving convection–diffusion equations, a difficult task is to construct discrete analogs of convective members that approximate continuous ones with the most excellent accuracy. A linear combination of the central difference scheme [19] and the Upwind Leapfrog scheme [20] was constructed to solve this task with weight coefficients selected to minimize the approximation error. The geometry of the computational domain is recorded using the method of considering the filling of rectangular cells with a material medium (this can be liquid, soil, and air), in particular, liquid. This makes it possible to increase the smoothness and accuracy of the finite–difference solution of hydrodynamic problems with a complex shape on the boundary surface. The central difference scheme for the transport equation considering the cell filling function [21] is written in the form:
( q 0 ) i , j c i , j n + 1 c i , j n τ + ( q 1 ) i , j u i + 1 / 2 , j c i + 1 , j n c i , j n 2 h x + ( q 2 ) i , j u i 1 / 2 , j c i , j n c i 1 , j n 2 h x = 0 .
The Upwind Leapfrog scheme for the transport equation is written as:
c i n + 1 c i n 2 τ + c i 1 n c i 1 n 1 2 τ + u i , j c i n c i 1 n h = 0 .
A scheme for the transfer equation taking into account the cell filling function, which is a linear combination of the central difference scheme and the Upwind Leapfrog scheme with weight coefficients 1 / 2 and ( q 2 ) i , j , can be written accordingly:
( q 0 ) i , j + ( q 2 ) i , j 2 c i , j n + 1 c i , j n τ + ( q 2 ) i , j c i 1 n c i 1 n 1 2 τ + ( q 1 ) i , j u i + 1 / 2 , j c i + 1 , j n c i , j n 4 h x + 5 ( q 2 ) i , j u i 1 / 2 , j c i , j n c i 1 , j n 4 h x = 0 .
Approximate Equations (5) and (6) using scheme (7) and taking into account the function of filling the cells [21]:
  • Difference scheme for Equation (5) describing transport along the direction Ox:
    ( q 0 ) i , j + ( q 2 ) i , j 2 c i , j n + 1 / 2 c i , j n τ + ( q 2 ) i , j c i 1 , j n 1 / 2 c i 1 , j n 1 2 τ + 5 ( q 2 ) i , j u i 1 / 2 , j c i , j n c i 1 , j n 4 h x + min ( ( q 1 ) i , j , ( q 2 ) i , j ) u i + 1 / 2 , j c i + 1 , j n c i , j n 4 h x = 3 ( q 1 ) i , j μ i + 1 / 2 , j c i + 1 , j n c i , j n 2 h x 2 3 ( q 2 ) i , j μ i 1 / 2 , j c i , j n c i 1 , j n 2 h x 2 | ( q 1 ) i , j ( q 2 ) i , j | μ i , j α x c i , j n + β x h x ,       u i , j 0 ;
    ( q 0 ) i , j + ( q 1 ) i , j 2 c i , j n + 1 / 2 c i , j n τ + ( q 1 ) i , j c i + 1 , j n 1 / 2 c i + 1 , j n 1 2 τ + 5 ( q 1 ) i , j u i + 1 / 2 , j c i + 1 , j n c i , j n 4 h x + min ( ( q 1 ) i , j , ( q 2 ) i , j ) u i 1 / 2 , j c i , j n c i 1 , j n 4 h x = 3 ( q 1 ) i , j μ i + 1 / 2 , j c i + 1 , j n c i , j n 2 h x 2 3 ( q 2 ) i , j μ i 1 / 2 , j c i , j n c i 1 , j n 2 h x 2 | ( q 1 ) i , j ( q 2 ) i , j | μ i , j α x c i , j n + β x h x ,       u i , j < 0 ;
  • Difference scheme for Equation (6) describing transport along the direction Oy:
    ( q 0 ) i , j + ( q 4 ) i , j 2 c i , j n + 1 c i , j n + 1 / 2 τ + ( q 4 ) i , j c i , j 1 n c i , j 1 n 1 / 2 2 τ + 5 ( q 4 ) i , j v i , j 1 / 2 c i , j n + 1 / 2 c i , j 1 n + 1 / 2 4 h y + + min ( ( q 3 ) i , j , ( q 4 ) i , j ) v i , j + 1 / 2 c i , j + 1 n + 1 / 2 c i , j n + 1 / 2 4 h y = 3 ( q 3 ) i , j μ i , j + 1 / 2 c i , j + 1 n + 1 / 2 c i , j n + 1 / 2 2 h y 2 3 ( q 4 ) i , j μ i , j 1 / 2 c i , j n + 1 / 2 c i , j 1 n + 1 / 2 2 h y 2 | ( q 3 ) i , j ( q 4 ) i , j | μ i , j α y c i , j n + 1 / 2 + β y h y ,       v i , j 0 ;
    ( q 0 ) i , j + ( q 3 ) i , j 2 c i , j n + 1 c i , j n + 1 / 2 τ + ( q 3 ) i , j c i , j + 1 n c i , j + 1 n 1 / 2 2 τ + 5 ( q 3 ) i , j v i , j + 1 / 2 c i , j + 1 n + 1 / 2 c i , j n + 1 / 2 4 h y + min ( ( q 3 ) i , j , ( q 4 ) i , j ) v i , j 1 / 2 c i , j n + 1 / 2 c i , j 1 n + 1 / 2 4 h y = 3 ( q 3 ) i , j μ i , j + 1 / 2 c i , j + 1 n + 1 / 2 c i , j n + 1 / 2 2 h y 2 3 ( q 4 ) i , j μ i , j 1 / 2 c i , j n + 1 / 2 c i , j 1 n + 1 / 2 2 h y 2 | ( q 3 ) i , j ( q 4 ) i , j | μ i , j α y c i , j n + 1 / 2 + β y h y ,       v i , j < 0 .
To solve the system of difference equations obtained as a result of discretizing a continuous model of the dynamics of phytoplankton populations (1)–(3), an adaptive modified alternating triangular iterative method was chosen. It requires the least number of iterations for a given accuracy when solving a skew-symmetric operator under the condition’s boundedness of the grid Péclet number [22].

2.3. Research of the Stability of a Central Difference Scheme and a Upwind Leapfrog Scheme Linear Combination

Write a scheme that is a linear combination of the central difference scheme and the Upwind Leapfrog scheme:
q i n + 1 q i n τ + q i 1 n q i 1 n 1 2 τ + u q i + 1 n + 4 q i n 5 q i 1 n 4 h = 0 .
Research the stability of (8) scheme by the harmonic method or the Neumann method [23]. Let q i n = φ n e j k i , where j = 1 , substitute in (8):
φ 1 τ + e j k e j k φ 2 τ + u e j k + 4 5 e j k 4 h = 0   or
φ 2 + ( ( 1 2 + u τ h ) ( 1 cos k ) 1 2 + j sin k 2 ( 1 + 3 u τ h sin k ) ) φ e j k 2 = 0 .
Solve the quadratic equation for φ :
φ 1 , 2 = ( 1 4 u τ h ) ( 1 cos k ) + 1 4 j sin k 4 ( 1 + 3 u τ h ) ± ( ( 1 4 u τ h ) ( 1 cos k ) + 1 4 j sin k 4 ( 1 + 3 u τ h ) ) 2 + e j k 2 .
The root φ 2 is not a solution, because φ 2 1 at τ = 0 . Therefore, let us move on to the replacement u τ / h = x and denote the absolute value of the function φ 1 ( x , k ) as:
ψ = | ( 1 4 x ) ( 1 cos k ) + 1 4 j sin k 4 ( 3 x 1 ) | + ( ( 1 4 x ) ( 1 cos k ) + 1 4 j sin k 4 ( 3 x 1 ) ) 2 + e j k 2 .
Study the behavior of the values of functions ψ ( x , k ) . Let us take meanings k [ 0 , 2 π ] and meanings x [ 0 , 1 ] . Figure 1 demonstrate that the values of ψ ( x , k ) belong to the segment [ 0 , 1 ] at k = π and x [ 0 , 0.75 ] . Under these conditions, the new scheme is stable.

2.4. Research of the Accuracy of a Central Difference Scheme and an Upwind Leapfrog Scheme Linear Combination

Consider the initial-boundary value problem for an equation of parabolic type with the lowest derivative:
q t + u q x = μ q x x ,
where t [ 0 , T ] , x [ 0 , L ] , u = c o n s t , μ = c o n s t , with initial conditions and boundary conditions q ( 0 , x ) = q 0 ( x ) , q ( t , 0 ) = q ( t , L ) = 0 .
Let us find an analytical solution to the problem (9) q ( x , t ) C 2 ( 0 < x < L ) C ( 0 x L ) C 1 ( 0 < t < + ) C ( 0 t < + ) .
Write the finite sum of the trigonometric Fourier series for the function q ( x , t ) in complex form:
q ( x , t ) = m = N N C m ( t ) exp ( j ω m x ) ,
where ω = π / L , m is the harmonic number; C m ( t ) = 2 L 0 L q ( x , t ) exp ( j ω m x ) d x is the complex amplitude of the m–th harmonic, j = 1 . Using a segment of a Fourier series (finite), we obtain an exact solution in the one-dimensional case. A finite trigonometric basis is used since the number of grid nodes is limited.
Substituting (10) into (9), obtains:
( m = N N C m ( t ) exp ( j ω m x ) ) t + u ( m = N N C m ( t ) exp ( j ω m x ) ) x = μ ( m = N N C m ( t ) exp ( j ω m x ) ) x x .
Change the sequence of operations of differentiation and summation of the series, calculate the derivative with respect to space:
m = N N ( C m ( t ) ) t exp ( j ω m x ) + j u ω m m = N N C m ( t ) exp ( j ω m x ) = μ ω 2 m 2 m = N N C m ( t ) exp ( j ω m x ) .
Considering that the functions exp(jωmx) are linearly independent, obtaining:
( C m ( t ) ) t = ( j u ω m + μ ω 2 m 2 ) C m ( t )
The solution to Equation (11) has the form:
C m ( t ) = C m ( 0 ) exp ( ( μ ω 2 m 2 + j u ω m ) t ) .
Substituting into (10) and taking into account the initial and boundary conditions, we obtain a solution to the Equation (9):
q ( x , t ) = m = N N C m ( 0 ) exp ( ( μ ω 2 m 2 + j u ω m ) t ) exp ( j ω m x ) .
Study the accuracy of the difference scheme considers a uniform space-time grid ω = ω τ × ω ¯ x , where ω ¯ x = { x i | x i = i h ;       i = 0 , 1 , , N x ;       N x h = L } , ω τ = { t n |   t n = n τ ;     n = 0 , 1 , , N t ;     N t τ = T } , τ is the time step, h is the space step, Nt is the number of time steps, Nx is the number of space steps.
Equation (10) is approximated as:
q i n + 1 q i n τ + q i 1 n q i 1 n 1 2 τ + u q i + 1 n + 4 q i n 5 q i 1 n 4 h 3 μ q i + 1 n 2 q i n + q i 1 n 2 h 2 = 0 ,
where:
q i = m = N N C m exp ( j m ω x i ) .
Substituting (14) into (13) and obtaining:
m = N N C m n + 1 exp ( j ω m h i ) C m n exp ( j ω m h i ) τ + m = N N C m n exp ( j m ω h ( i 1 ) ) C m n 1 exp ( j m ω h ( i 1 ) ) 2 τ + + u m = N N C m n exp ( j ω m h ( i + 1 ) ) + 4 C m n exp ( j ω m h i ) 5 C m n exp ( j ω m h ( i 1 ) ) 4 h 3 μ m = N N C m n exp ( j ω m h ( i + 1 ) ) 2 C m n exp ( j ω m h i ) + C m n exp ( j ω m h ( i 1 ) ) 2 h 2 = 0 .
Doing simple arithmetic operations, obtaining:
m = N N C m n + 1 C m n τ exp ( j ω m h i ) + m = N N C m n C m n 1 2 τ exp ( j m ω h ) exp ( j m ω h i ) + u m = N N C m n exp ( j ω m h ) + 4 C m n 5 C m n exp ( j ω m h ) 4 h exp ( j ω m h i ) 3 μ m = N N C m n exp ( j ω m h ) 2 C m n + C m n exp ( j ω m h ) 2 h 2 exp ( j ω m h i ) = 0 .
Due to linear independence exp ( j m ω i ) write:
C m n + 1 C m n τ + C m n C m n 1 2 τ exp ( j ω m h ) + u C m n exp ( j ω m h ) + 4 5 exp ( j ω m h ) 4 h 3 μ C m n exp ( j ω m h ) 2 + exp ( j ω m h ) 2 h 2 = 0 .
For τ 0 from (15) it follows:
( C m ( t ) ) t ( 1 + exp ( j ω m h ) 2 ) + u C m n exp ( j ω m h ) + 4 5 exp ( j ω m h ) 4 h 3 μ C m n exp ( j ω m h ) 2 + exp ( j ω m h ) 2 h 2 = 0 ;
( C m ( t ) ) t = ( u exp ( j ω m h ) + 4 5 exp ( j ω m h ) 2 h ( 2 + exp ( j ω m h ) ) + μ 3 ( exp ( j ω m h ) 2 + exp ( j ω m h ) ) h 2 ( 2 + exp ( j ω m h ) ) ) C m n .
It can be seen from the last expression that the continuous problem (9) differs from the discrete one (13) by the quantities α 1 = 1 exp ( j ω m h ) + 4 5 exp ( j ω m h ) 2 j ω m h ( 2 + exp ( j ω m h ) ) and α 2 = 1 3 ( exp ( j ω m h ) + 2 exp ( j ω m h ) ) ω 2 m 2 h 2 ( 2 + exp ( j ω m h ) ) , accordingly.
When the order of the error in approximating the convective and diffusion terms in space is studied, we find that scheme (13) approximates the convective term with the third order of accuracy in space, and the diffusion term with the first one.
Comment. The quantity r = π / ω m h describes the number of nodes on half of the wave period while the estimate π > ω m h takes place. Hence it follows that the accuracy of the solution depends on the number of nodes on half of the wave period.
Figure 2 show a graph of the function α 1 ( r ) , describing the dependence of the approximation error of the convective and diffusion operators, respectively, by the difference scheme (13) on the number of nodes used to solve the problem (9) in comparison with the central difference scheme.
Analyzing the graphs, we can conclude that scheme (13), a linear combination of the central difference scheme and the Upwind Leapfrog scheme, effectively solves convection–diffusion problems in which convective transfer prevails over diffusion and the values of the grid Péclet number are in the range 2 < P e 20 .

2.5. Software Implementation

For the mathematical modeling of the phytoplankton development dynamics, taking into account the transformation of forms of phosphorus, nitrogen, and silicon, a software package (SP) has been developed. The SP simulates the dynamics of the development of three main types of summer phytoplankton—blue-green algae (Aphanizomenon flos-aquae), green algae (Chlorella vulgaris), diatoms (Skeletonema costatum), their competition for nutrients and the transformation of the forms of these nutrients—phosphorus, nitrogen, and silicon, their absorption, excretion, and transition from one biochemical compound to another, the influence of salinity, and temperature on the growth rate of phytoplankton, are all taken into account. The developed software makes it possible to predict the dynamics of the ecosystem’s development of shallow water bodies using the example of the Azov Sea and the Taganrog Bay with an increase in water salinity.
The program block “Biogeochemical Cycles” uses a three-dimensional vector of the speed of movement of the water flow in the Azov Sea, as well as three-dimensional fields of salinity and temperature, calculated using the program module “Calculation of the aquatic environment movement” (“Azov3D.exe”). In the software module, the calculation took into account such factors and processes as the Coriolis force, turbulent exchange, complex geometry of the bottom and coastline, evaporation, river flows, wind currents, and friction on the bottom [24]. The scheme of the algorithm of the program module “Calculation of the aquatic environment movement” is shown in Figure 3.

3. Results

Modeling is performed in a rectangular area, the dimensions of which correspond to the physical dimensions of the Azov Sea, using a uniform grid. The time interval is 30 days, the values of the temperature field are taken in accordance with the long-term average data for July. Initial concentration values of blue-green algae (Aphanizomenon flos-aquae)—2.6 mg/L, green algae (Chlorella vulgaris)—2.5 mg/L, diatoms (Skeletonema costatum)—are 0.9 mg/L, distributions are uniform. The optimal salinity for the first two species of phytoplankton—6‰, for the third—are 12‰, the coefficients of the width of the salinity tolerance interval are b 1 , 2 , 3 = 2 . Figure 4 show the graphs of the dependencies of the growth rates of three phytoplankton species on salinity and the salinity field reconstructed from cartographic information based on high-order approximation schemes [25] for the normal salinity level.
A series of experiments were carried out (Figure 5 and Figure 6), in each of which the salinity level increased. The initial salinity level was taken corresponding to the data of long-term observations in full-flowing years, when a sufficient amount of freshwater, including the Don and Kuban rivers, enters the water area of the Azov Sea.
According to literary sources and expeditionary studies, the normal salinity level of the Azov Sea corresponds to the following values: 0–9‰ for the Taganrog Bay, 11–12‰ in the main water area, and 12–14‰ closer to the Kerch Strait [26]. Figure 5 and Figure 6 show the results of modeling the dynamics of blue-green and green algae, respectively, in the Azov Sea and the Taganrog Bay at different levels of water salinity; isohaline (isolines connecting places of a reservoir with equal salinity) of salinity are indicated by a dotted line. In the experiments, the salinity at the outlet of the bay (a conventional line connecting the Dolzhanskaya and Belosaraiskaya spits) varies from 9‰ to 13.5‰.
Numerical experiments have shown that with an increase in salinity in the Azov Sea, blue-green and green algae move to the mouth of the Don River and their concentration decreases. This entails a change in the food pyramid base and the habitat of freshwater aquatic organisms of higher trophic levels—zooplankton, molluscs, and commercial fish species.

4. Discussion

In this paper, a three-dimensional mathematical model of the dynamics of phytoplankton populations is proposed, combined with a three-dimensional model of the hydrodynamics of shallow water bodies. The Azov3D software package implements a mathematical model of hydrodynamic processes in coastal systems based on the equations of motion (Navier–Stokes) in three coordinate directions and the continuity equation for a fluid with variable density. In contrast to the known mathematical models used in the works of other researchers, built on the Saint-Venant equations [27], this approach allows the consideration of many variables, such as the transport of heat and salts, the variable density of the aqueous medium, friction on the bottom and wind stresses, turbulent exchange, the Coriolis force, river runoff, evaporation, and precipitation. In addition, it is also considered the complicated and dynamically changing geometry of the waterbody, which makes it possible to increase the accuracy of modeling, detect different-scale eddy structures of currents in coastal systems, and act as natural traps for pollution.
The developed software package, built on the basis of models of hydrodynamics and biogeochemical cycles, showed stability at depth differences of 40–50 times, which is typical for coastal systems. Similarly, the known models (MARS 3D [9,10,11], POM, etc.) turned out to be computationally unstable. In contrast to the known models and programs, the developed software package made it possible to detect eddy structures (the Azov and Mediterranean seas), zones of hypoxia, and anaerobic contamination. In addition, it can reconstruct the ecological catastrophe in the Azov Sea and predict with high accuracy the forecast of an extreme storm that happened on 23–24 September 2014 in the Taganrog Bay, when, with a wind speed of 42 m/s and an average depth of 2 m, the level rise was more than 4 m.
The three-dimensional mathematical model of biogeochemical cycles, combined with the model of hydrodynamics, used in this work, makes it possible to predict the spatial-temporal dynamics of the main phytoplankton populations more accurately compared to the two-dimensional models discussed above [12,13,15,16], taking into account the changing hydrological regime, their consumption of nutrients, the transition of nutrients from one form to another. In contrast to the two-dimensional model used in [14], the three-dimensional model of hydrodynamics makes it possible to more accurately model wave processes under conditions of a complex, dynamically changing the geometry of the reservoir, friction on the bottom and wind stress on the free surface, turbulent and advective heat and mass transfer in three coordinate directions, Coriolis forces, river flows, precipitation and evaporation, heat and salt transportation.
The application of classical center-difference schemes for the approximation of advection–diffusion reaction problems, in the case of large values of the grid Péclet number, leads to the instability of the solution. When using the CABARET-type (Upwind Leapfrog) schemes described in work [28], oscillations occur. In this paper, a linear combination of these difference schemes with weight coefficients was obtained to minimize the approximation error. This difference scheme, while solving diffusion–convection problems, does not have grid viscosity. Consequently, this scheme more accurately describes the behavior of the solution in the case of large grid Péclet numbers (2 < Peh ≤ 20). It also preserves the smoothness of the solution at the interface between the media boundaries when solving dynamic problems with a complex shape of the boundary surface (there are no oscillations associated with the stepwise approximation of the borders). The study of the proposed linear combination of difference schemes has shown that this scheme allows one to approximate the convective term with the third order of accuracy. These schemes in solving hydrodynamic problems make it possible to describe small-sized eddys that arise in the coastal part of shallow water bodies more accurately.
The developed software is built on interconnected models that more accurately describe biogeochemical processes, the distribution of nutrient concentrations in shallow water bodies under conditions of significantly changing salinity, and the density of the aquatic environment compared with the known models. Furthermore, considering several hydrodynamic and hydrophysical factors, the proposed model makes it possible to increase the accuracy of modeling the transport processes of nutrients in coastal systems, such as the Azov Sea, in complex spatially heterogeneous structures of currents, a complex coastline, and bottom topography. The simulation data are consistent with space sensing data and the available experimental data from previous periods [29].

5. Conclusions

As a result of this research, the influence of salinity as an abiotic factor concerning the dynamics of the development of the main species of phytoplankton populations in the Azov Sea and Taganrog Bay, was studied. A numerical experiment carried out on the basis of a scenario approach helped to track the geographical dynamics of green and blue-green algae and the displacement of their habitats to the mouth of the Don River with an increase in salinity due to a decrease in freshwater river runoff. This results in a change in the habitats of aquatic organisms of higher trophic levels—zooplankton, commercial, and valuable fish species, in which freshwater species are replaced by marine species. The use of a linear combination of central and Upwind Leapfrog difference schemes for discretization made it possible to increase the accuracy of modeling at large values of the grid Péclet number (Peh > 2). The developed software package includes modules for calculating: a three-dimensional vector of the speed of movement of a water flow in the Azov Sea, three-dimensional fields of salinity and temperature, and concentrations of the main species of summer phytoplankton at different levels of sea salinity. The developed software also makes it possible to predict the geographical dynamics of the main phytoplankton species with an increase in the salinity of the Azov Sea and the Taganrog Bay due to a decrease in the fresh runoff of the Don and Kuban rivers and the inflow of saline waters of the Black Sea with a dynamically changing wind regime through the Kerch Strait.

Author Contributions

Conceptualization, A.S. and A.C.; methodology, A.C.; software, Y.B.; validation, A.S. and Y.B.; formal analysis, A.B.; investigation, A.C. and Y.B.; resources, A.S. and Y.B.; data curation, A.S.; writing—original draft preparation, A.C. and Y.B.; writing—review and editing, A.S. and A.B.; visualization, Y.B.; supervision, A.S.; project administration, A.B. and B.M.; funding acquisition, A.B. and B.M. All authors have read and agreed to the published version of the manuscript.

Funding

The study was supported by a grant from the Russian Science Foundation (project No. 21-71-20050).

Institutional Review Board Statement

The study does not include humans or animal research.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Acknowledgments

The authors would like to acknowledge the administration of Don State Technical University for resource and financial support.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Matishov, G.G.; Kovaleva, G.V.; Yasakova, O.N. Abnormal salinization in the Taganrog estuary and the Don delta. Sci. South Russ. 2016, 12, 43–50. [Google Scholar]
  2. Falkowski, P.G.; Oliver, M.J. Mix and match: How climate selects phytoplankton. Nat. Rev. Microbiol. 2007, 5, 813–819. [Google Scholar] [CrossRef] [PubMed]
  3. Ilyichev, V.G. Ideas of evolutionary ecology in models of aquatic ecological systems. Water Resour. 1993, 29, 5–11. [Google Scholar]
  4. Menshutkin, V.V.; Filatov, N.N.; Potakhin, M.S. Expert system “Lakes of Karelia”. 2. Classification of lakes. Water Resour. 2009, 36, 300–311. [Google Scholar]
  5. O’Reilly, J.E.; Zetlin, C. Seasonal, Horizontal, and Vertical Distribution of Phytoplankton Chlorophyll a in the Northeast U.S. Continental Shelf Ecosystem; U.S. Department of Commerce: Seattle, WA, USA, 1998; p. 119.
  6. Abrosov, N.S.; Bogolyubov, A.A. Ecological and Genetic Patterns of Coexistence and Co-Evolution of Species; Nauka: Novosibirsk, Russia, 1988; p. 327. [Google Scholar]
  7. Kovaleva, G.V. Check-list of benthic and plankton microalgae of a coastal part of the Azov Sea and adjoining reservoirs. In Proceedings of the International Scientific Conference and VIIth Marine Biology School, Rostov-on-Don, Russia, 9–13 June 2008; Matishov, G., Ed.; SSC RAS Publishers: Rostov-on-Don, Russia, 2008. [Google Scholar]
  8. Vinberg, G.G. Energy principle of studying trophic connections and productivity of ecological systems. Zool. J. 1962, 41, 1618–1630. [Google Scholar]
  9. MARS. Système de Modélisation de l’Environnement Côtier. Available online: https://wwz.ifremer.fr/mars3d (accessed on 10 September 2021).
  10. Alekseenko, E.; Roux, B. Contribution to remediation of brackish lagoon: 3D simulation of salinity, bottom currents and resuspension of bottom sediments by strong winds. Estuar. Coast. Shelf Sci. 2019, 216, 27–37. [Google Scholar] [CrossRef]
  11. Ménesguen, A.; Dussauze, M.; Dumas, F.; Thouvenin, B.; Garnier, V.; Lecornu, F.; Répécaud, M. Ecological model of the Bay of Biscay and English Channel shelf for environmental status assessment part 1: Nutrients, phytoplankton and oxygen. Ocean Model. 2019, 133, 56–78. [Google Scholar] [CrossRef] [Green Version]
  12. Weng, X.; Jiang, C.; Zhang, M.; Yuan, M.; Zeng, T. Numeric Study on the Influence of Sluice-Gate Operation on Salinity, Nutrients and Organisms in the Jiaojiang River Estuary, China. Water 2020, 12, 2026. [Google Scholar] [CrossRef]
  13. Zhu, C.; Liang, Q.; Yan, F.; Hao, W. Reduction of Waste Water in Erhai Lake Based on MIKE21 Hydrodynamic and Water Quality Model. Sci. World J. 2013, 2013, 958506. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Sohrt, V.; Hein, S.S.V.; Nehlsen, E.; Strotmann, T.; Fröhle, P. Model Based Assessment of the Reflection Behavior of Tidal Waves at Bathymetric Changes in Estuaries. Water 2021, 13, 489. [Google Scholar] [CrossRef]
  15. Yakushev, E.V.; Wallhead, P.; Renaud, P.E.; Ilinskaya, A.; Protsenko, E.; Yakubov, S.; Pakhomova, S.; Sweetman, A.K.; Dunlop, K.; Berezina, A.; et al. Understanding the biogeochemical impacts of fish farms using a benthic-pelagic model. Water 2020, 12, 2384. [Google Scholar] [CrossRef]
  16. Yakushev, E.V.; Makkaveev, P.N.; Polukhin, A.A.; Protsenko, E.A.; Stepanova, S.V.; Khlebopashev, P.V.; Yakubov, S.K.; Staalstrøm, A.; Norli, M. Hydrochemical studies in coastal waters of the Spitsbergen Archipelago in 2014–2015. Oceanology 2016, 56, 763–765. [Google Scholar] [CrossRef]
  17. Tyutyunov, Y.; Senina, I.; Arditi, R. Clustering due to acceleration in the response to population gradient: A simple self-organization model. Am. Nat. 2004, 164, 722–735. [Google Scholar] [CrossRef] [PubMed]
  18. Silkin, V.A.; Abakumov, A.I.; Pautova, L.A.; Pakhomova, S.V.; Lifanchuk, A.V. Mechanisms of regulation of invasive processes in phytoplankton on the example of the north-eastern part of the Black Sea. Aquat. Ecol. 2016, 50, 221–234. [Google Scholar] [CrossRef]
  19. Wilkes, M.V. A Short Introduction to Numerical Analysis; Cambridge University Press: Cambridge, UK, 1966. [Google Scholar]
  20. Thomas, J.P.; Roe, P.L. Development of non-dissipative numerical schemes for computational aeroacoustics. In Proceedings of the 11th AIAA Computational Fluid Dynamics Conference, Orlando, FL, USA, 6–9 July 1993. [Google Scholar]
  21. Sukhinov, A.I.; Chistyakov, A.E.; Protsenko, E.A.; Sidoryakina, V.V.; Protsenko, S.V. Accounting Method of Filling Cells for the Solution of Hydrodynamics Problems with a Complex Geometry of the Computational Domain. Math. Models Comput. Simul. 2020, 12, 232–245. [Google Scholar] [CrossRef]
  22. Sukhinov, A.I.; Chistyakov, A.E.; Belova, Y.V. The difference scheme for the two-dimensional convection-diffusion problem for large peclet numbers. MATEC Web Conf. 2018, 226, 04030. [Google Scholar] [CrossRef] [Green Version]
  23. Neumann, J.; Richtmyer, R.D. On the stability of finite difference matrices. J. Soc. Industr. Appl. Math. 1965, 2, 119–128. [Google Scholar]
  24. Sukhinov, A.I.; Chistyakov, A.E.; Alekseenko, E.V. Numerical realization of the three-dimensional model of hydrodynamics for shallow water basins on a high-performance system. Math. Models Comput. Simul. 2011, 3, 562–574. [Google Scholar] [CrossRef]
  25. Sukhinov, A.I.; Belova, Y.V.; Filina, A.A. Parallel implementation of substance transport problems for restoration the salinity field based on schemes of high order of accuracy. In Proceedings of the CEUR Workshop Proceedings, Otzenhausen, Germany, 8–12 September 2019; p. 2500. [Google Scholar]
  26. Matishov, G.G. Ecological Atlas of the Azov Sea; Publishing House of the SSC RAS: Rostov-on-Don, Russia, 2011; p. 328. [Google Scholar]
  27. Chikin, A.L.; Kleshchenkov, A.V.; Chikina, L.G. Simulating salinity variations in the Gulf of Taganrog at storm surges. Water Resour. 2019, 46, 919–925. [Google Scholar] [CrossRef]
  28. Glotov, V.Y.; Goloviznin, V.M. CABARET scheme in velocity-pressure formulation for two-dimensional incompressible fluids. Comput. Math. Math. Phys. 2013, 53, 721–735. [Google Scholar] [CrossRef]
  29. Terziev, F.S. (Ed.) Hydrometeorology and Hydrochemistry of the USSR Seas. In The Azov Sea; Gidrometeoizdat: St. Petersburg, Russia, 1991; Volume 5, p. 237. [Google Scholar]
Figure 1. The function ψ ( x , π ) graph.
Figure 1. The function ψ ( x , π ) graph.
Mathematics 09 03025 g001
Figure 2. The graph of the dependence of the convective operator approximation error on the number of nodes: 1—for the scheme (13); 2—for the central difference scheme.
Figure 2. The graph of the dependence of the convective operator approximation error on the number of nodes: 1—for the scheme (13); 2—for the central difference scheme.
Mathematics 09 03025 g002
Figure 3. Program module block diagram.
Figure 3. Program module block diagram.
Mathematics 09 03025 g003
Figure 4. The salinity influence on the phytoplankton growth: (a) graphs of the growth rate of 1—green, 2—blue-green, 3—diatoms on salinity; (b) reconstructed salinity field.
Figure 4. The salinity influence on the phytoplankton growth: (a) graphs of the growth rate of 1—green, 2—blue-green, 3—diatoms on salinity; (b) reconstructed salinity field.
Mathematics 09 03025 g004
Figure 5. Distributions of blue-green algae (Aphanizomenon flos-aquae) concentrations at different salinity levels: (a) initial; (b) increased by 10%; (c) increased by 30%; (d) increased by 50%.
Figure 5. Distributions of blue-green algae (Aphanizomenon flos-aquae) concentrations at different salinity levels: (a) initial; (b) increased by 10%; (c) increased by 30%; (d) increased by 50%.
Mathematics 09 03025 g005aMathematics 09 03025 g005b
Figure 6. Distributions of green algae (Chlorella vulgaris) concentrations at different salinity levels: (a) initial; (b) increased by 10%; (c) increased by 30%; (d) increased by 50%.
Figure 6. Distributions of green algae (Chlorella vulgaris) concentrations at different salinity levels: (a) initial; (b) increased by 10%; (c) increased by 30%; (d) increased by 50%.
Mathematics 09 03025 g006
Table 1. Nutrients in the model dynamics of phytoplankton.
Table 1. Nutrients in the model dynamics of phytoplankton.
NumberDesignationTitle
1F1green algae Chlorella vulgaris
2F2blue-green algae Aphanizomenon flos-aquae
3F3diatom algae Skeletonema costatum
4PO4phosphates
5POPsuspended organic phosphorus
6DOPdissolved organic phosphorus
7NO3nitrates
8NO2nitrites
9NH4ammonium
10Sidissolved inorganic silicon (silicic acids)
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sukhinov, A.; Belova, Y.; Chistyakov, A.; Beskopylny, A.; Meskhi, B. Mathematical Modeling of the Phytoplankton Populations Geographic Dynamics for Possible Scenarios of Changes in the Azov Sea Hydrological Regime. Mathematics 2021, 9, 3025. https://doi.org/10.3390/math9233025

AMA Style

Sukhinov A, Belova Y, Chistyakov A, Beskopylny A, Meskhi B. Mathematical Modeling of the Phytoplankton Populations Geographic Dynamics for Possible Scenarios of Changes in the Azov Sea Hydrological Regime. Mathematics. 2021; 9(23):3025. https://doi.org/10.3390/math9233025

Chicago/Turabian Style

Sukhinov, Alexander, Yulia Belova, Alexander Chistyakov, Alexey Beskopylny, and Besarion Meskhi. 2021. "Mathematical Modeling of the Phytoplankton Populations Geographic Dynamics for Possible Scenarios of Changes in the Azov Sea Hydrological Regime" Mathematics 9, no. 23: 3025. https://doi.org/10.3390/math9233025

APA Style

Sukhinov, A., Belova, Y., Chistyakov, A., Beskopylny, A., & Meskhi, B. (2021). Mathematical Modeling of the Phytoplankton Populations Geographic Dynamics for Possible Scenarios of Changes in the Azov Sea Hydrological Regime. Mathematics, 9(23), 3025. https://doi.org/10.3390/math9233025

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop