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

: Increased inﬂuence 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 inﬂuence 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 (Pe h > 2). A software package has been developed that implements interrelated models of hydrodynamics and biogeochemical cycles. A modiﬁed 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.


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

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.

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]: The following equations describe chemical-biological sources: 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 PD is the specific speed of autolysis POP; K PN is the phosphatification coefficient POP; K DN 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 Si are the normalization coefficients between the content of N, P, Si in organic matter. The following expressions determine the growth rate of phytoplankton: where K NF is the maximum specific growth rate of phytoplankton.
Functions of the dependence of the growth rate of aquatic organisms on temperature and salinity: where k s = 1, T opt , S opt 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 opt and S opt 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: - where K NO 3 is the half-saturation constant of nitrates, K NH 4 is the half-saturation constant of ammonium, K psi 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 u n 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 c i : where ε i are the non-negative constants, i ∈ M; ε i 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 c 0 i (x, y, z) are the preset functions, fields of water flow vector velocity, salinity, and temperature are input data for this model.

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, 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): 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 1/2 and (q 2 ) i,j , 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].

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 q n i = ϕ n · e jki , where j = √ −1, substitute in (8): Solve the quadratic equation for ϕ: 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: 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.

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 Let us find an analytical solution to the problem (9) , Write the finite sum of the trigonometric Fourier series for the function ( ) , q x t in complex form:

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 t ∈ [0, T], x ∈ [0, L], u = const, µ = const, 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 < Write the finite sum of the trigonometric Fourier series for the function q(x, t) in complex form: where ω = π/L, m is the harmonic number; C m (t) = 2 L L 0 q(x, t) exp(jωmx)dx 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: 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- . . , N t ; N t τ = T}, τ is the time step, h is the space step, N t is the number of time steps, N x is the number of space steps.
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ωmh)+4−5 exp(−jωmh) 2jωmh(2+exp(−jωmh)) 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 r = π/ωmh describes the number of nodes on half of the wave period while the estimate π > ωmh 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 convectiondiffusion problems in which convective transfer prevails over diffusion and the values of the grid Péclet number are in the range 2 < Pe ≤ 20.

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 discrete one (13) by the quantities 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 mh π ω = describes the number of nodes on half of the wave period while the estimate mh π ω > 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 2 0 Pe < ≤ . 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.

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.

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

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 ( Figures 5 and 6), in each of which the salinity level increased. The initial salinity level was taken corresponding to the data of longterm 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. A series of experiments were carried out ( Figures 5 and 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. A series of experiments were carried out ( Figures 5 and 6), in each of which the salinity level increased. The initial salinity level was taken corresponding to the data of longterm 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]. Figures 5 and 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. 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]. Figures 5 and 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. 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]. Figures 5 and 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 freshwa-ter aquatic organisms of higher trophic levels-zooplankton, molluscs, and commercial fish species.

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-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 spatialtemporal 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 advectiondiffusion 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 < Pe h ≤ 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].

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 (Pe h > 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, threedimensional 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.