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/m
3. 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 (NH
4, NO
3, PO
4, 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]:
where
is the concentration of
i-th component, [mg/L];
i ∈
M,
M = {F
1, F
2, F
3, PO
4, POP, DOP, NO
3, NO
2, NH
4, Si};
are the components of water flow velocity vector, [m/s];
k is the turbulent exchange coefficient, [m
2/s];
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
indicates the type of substance (
Table 1).
The following equations describe chemical-biological sources:
where
is the specific respiration rate of phytoplankton;
is the specific rate of phytoplankton dying;
is the specific rate of phytoplankton excretion;
is the specific speed of autolysis POP;
is the phosphatification coefficient POP;
is the phosphatification coefficient DOP;
is the specific rate of oxidation of ammonium to nitrites in the process of nitrification;
is the specific rate of oxidation of nitrites to nitrates in the process of nitrification;
,
,
are the normalization coefficients between the content of N, P, Si in organic matter.
The following expressions determine the growth rate of phytoplankton:
where
is the maximum specific growth rate of phytoplankton.
Functions of the dependence of the growth rate of aquatic organisms on temperature and salinity:
where
,
,
are the temperature and salinity optimal for a given type of aquatic organisms; and
,
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.
and
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 ,
where is the half-saturation constant of phosphates.
- -
for silicon ,
where is the half-saturation constant of silicon.
- -
for nitrogen
,
where
is the half-saturation constant of nitrates,
is the half-saturation constant of ammonium,
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:
where
are the non-negative constants,
i ∈
M;
take into account the sinking of algae to the bottom and their flooding for
i ∈ {F
1, F
2, F
3} and take into account the absorption of nutrients by bottom sediments for
i ∈ {PO
4, POP, DOP, NO
3, NO
2, NH
4, Si}.
It is also necessary to add the initial conditions:
where
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.
with boundary conditions:
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:
where
τ is the time step,
,
are the space steps,
is the upper time limit,
,
are the space boundaries,
,
are the characteristic dimensions of the computational grid.
The splitting schemes in space are used to discretize the homogeneous Equation (4):
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:
The Upwind Leapfrog scheme for the transport equation is written as:
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
and
, can be written accordingly:
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:
Difference scheme for Equation (6) describing transport along the direction
Oy:
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:
Research the stability of (8) scheme by the harmonic method or the Neumann method [
23]. Let
, where
, substitute in (8):
Solve the quadratic equation for
:
The root
is not a solution, because
at
. Therefore, let us move on to the replacement
and denote the absolute value of the function
as:
Study the behavior of the values of functions
. Let us take meanings
and meanings
.
Figure 1 demonstrate that the values of
belong to the segment
at
and
. 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:
where
,
,
,
, with initial conditions and boundary conditions
,
.
Let us find an analytical solution to the problem (9) .
Write the finite sum of the trigonometric Fourier series for the function
in complex form:
where
,
m is the harmonic number;
is the complex amplitude of the
m–th harmonic,
. 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:
Change the sequence of operations of differentiation and summation of the series, calculate the derivative with respect to space:
Considering that the functions exp(
jωmx) are linearly independent, obtaining:
The solution to Equation (11) has the form:
Substituting into (10) and taking into account the initial and boundary conditions, we obtain a solution to the Equation (9):
Study the accuracy of the difference scheme considers a uniform space-time grid , where , , τ 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:
where:
Substituting (14) into (13) and obtaining:
Doing simple arithmetic operations, obtaining:
Due to linear independence
write:
For
from (15) it follows:
It can be seen from the last expression that the continuous problem (9) differs from the discrete one (13) by the quantities and , 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 describes the number of nodes on half of the wave period while the estimate 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
, 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.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
.
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].