Sufﬁcient Conditions for the Existence and Uniqueness of the Solution of the Dynamics of Biogeochemical Cycles in Coastal Systems Problem

: The article considers a three-dimensional mathematical model of population dynamics based on a system of non-stationary parabolic advection-diffusion-reaction equations with lower derivatives describing the advective motion of the aquatic environment and non-linear source functions. In contrast to the previous authors’ works devoted to the description of this model and its numerical implementation, this article presents the results of an analytical study of the initial-boundary value problem corresponding to this model. For these purposes, the original initial-boundary value problem is linearized on a single time grid—for all nonlinear sources, their ﬁnal spatial distributions for the previous time step are used. As a result, a chain of initial-boundary value problems is obtained, connected by initial—ﬁnal data at each step of the time grid. For this chain of linearized problems, the existence and uniqueness of the solution of the initial-boundary value problem for the system of partial differential equations in the Hilbert space were researched. Numerical experiments were performed for model problems of the main types of phytoplankton populations in coastal systems under the inﬂuence of dynamically changing biotic and abiotic factors, the results of which are consistent with real physical experiments. The developed model, including the proposed model of biological kinetics, allows for the study of the productive and destructive processes of the shallow water body biocenosis to assess the state of the processes of reproduction of valuable and commercial ﬁsh participating in the food chain with selected species of summer phytoplankton.


Introduction
In recent decades, there has been a negative trend towards eutrophication of the South of Russia marine system waters, causing a rapid growth in phytoplankton populations, many of which are harmful and toxic. This, in turn, leads to a change in the species' composition and geography of the location of plankton populations, which are the basis of the food pyramid of coastal systems, such as the Taganrog Bay and the Azov Sea, the Tsimlyansk Reservoir, etc., degradation of separate components of the ecosystem, as well as entire communities of organisms in them. Prediction of such evolution of water systems, primarily associated with significant fluctuations in freshwater runoff, involves the further development of mathematical model methods and effective numerical methods for their implementation, which allow "losing" various scenarios for the development of the ecosystem of a water body based on non-stationary spatially-heterogeneous interconnected models of geochemical cycles and biological kinetics. order finite difference numerical method is used for the reformulated planetary geostrophic equations with an inviscid balance equation in article [25].
Earlier, in the article [26], the team of authors presented a description of some aspects of the construction and numerical implementation of the model of biogeochemical cycles and the biological kinetics of the multispecies populations model. However, the issues of analytical study of the model, related to the correctness of its formulation, remained in the shadows. This article, in which, along with the study of the correctness of the formulation of the initial-boundary value problem, computational experiments with an improved model are performed, eliminating this gap.
The considering model is based on a system of parabolic type equations with lower derivatives and non-linear source functions, which takes into account such important characteristics as river flows and sea currents, microturbulent vertical diffusion, flooding of matter and gravitational sedimentation, the spatially uneven distribution of temperature and salinity, as well as the interaction of the main biogenic substances-compounds of nitrogen, phosphorus, and the main types of plankton populations, including their growth, reproduction, natural decrease in abundance, etc. For the proposed model, nonlinear source functions are linearized on a uniform time grid, when the values of the nonlinear terms are determined as their final values on the previous time layer (with a delay), and a chain of initial and final solutions of the initial-boundary diffusion-convection-reaction problems is formed. The purpose of this research is to determine the sufficient conditions for the existence and uniqueness of solutions to the initial-boundary value problems of the planktonic populations and biogeochemical cycles dynamics. To do this, quadratic functionals are constructed, and inequalities are determined that guarantee their positivity and the existence and uniqueness of solutions using the energy method and the Gauss and Poincaré theorems. The research scientific novelty is the analytical study of the geochemical cycles and biological kinetics model, the linearization of the continuous model, the determination of the conditions for the positive definiteness of the operator of the equations system in the Hilbert space. This, in turn, guarantees the uniqueness of the problem solution and the unique and continuous dependence on the right-side function.
The Azov Sea, a unique water body in the South of Russia, has been chosen as a modeling object, for which this model was verified [17,[27][28][29][30]. The choice of this water body for the study is not accidental. The Azov Sea has a set of unique features, such as large fresh river water runoff, a large sea water salinity gradient, an abundance of nutrients coming with river runoff, temperature fluctuations greatly affecting the state of the sea ecosystem due to its shallow depth. These and other factors determine the biological diversity and high primary productivity of the reservoir. However, in recent decades, the frequency and severity of the consequences of adverse and dangerous phenomena in the Azov Sea were associated with eutrophication of the reservoir, abundant flowering of poisonous algae, the influx of pollutants, and the formation of extensive zones of hypoxia and anaerobicity. Almost every year, zones of deficiency of dissolved oxygen at the bottom are recorded in the Azov Sea. The inflow of pollutants with river runoff, the formation of a large amount of detritus due to abundant phytoplankton blooms in the warm season in the presence of different-scale whirlpools in the water flow, as a rule, lead to the appearance of hypoxia zones. The size and depth of these zones change annually, but with a significant and prolonged oxygen lack, so-called fish kill phenomena occur, when almost all benthic fauna, including fish, die, which has an extremely negative effect on the reproduction and fish productivity of the reservoir and entail the loss of aquatic biological resources.
In connection with the foregoing, the problem of constructing and researching integrated models of geochemical cycles and biological kinetics, including the conditions for the existence and uniqueness of solutions, is very relevant and has both fundamental scientific and applied significance.

Materials and Methods
A mathematical model of biological kinetics is used to study nonlinear effects in the dynamics of the most common phytoplankton species in the summer period.

Mathematical Model of the Phytoplankton Populations Dynamics
The model of the phytoplankton populations dynamics is based on a system of nonstationary equations of the convection-diffusion-reaction of parabolic type with nonlinear functions of sources and lower-order derivatives, which have the form: Chemical-biological sources are described by the following equations: 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 growth rate of phytoplankton is determined by the expressions: where K NF is the maximum specific growth rate of phytoplankton. Functions of the dependence of the growth rate of aquatic organisms on temperature (T) and salinity (S): where k s = 1; T opt , S opt are the temperature and salinity optimal for a given type of aquatic organisms; a l > 0, b l > 0 are the coefficients of the width of the range of aquatic organisms' tolerance to temperature and salinity, respectively. Functions describing biogen content: for phosphorus, where K PO 4 is the half-saturation constant of phosphates; -for silicon, where K Si is the half-saturation constant of silicon; -for nitrogen, 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. It is assumed that the coefficients used in the right-side functions are positive constants.
An initial-boundary value problem is posed in a cylindrical domain G for 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 be 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 q i : where ε i are the non-negative constants; i ∈ M; ε i consider the sinking of algae to the bottom and their flooding for i ∈ {F 1 , F 2 , F 3 } and consider the absorption of nutrients by bottom sediments for i ∈ {PO 4 , POP, DOP, NO 3 , NO 2 , NH 4 , Si}.
Add the initial values for the studied substances, as well as the water flow velocity vector, salinity, and temperature fields at any time for the system of Equation (1): T(x, y, z, 0) = T 0 (x, y, z), S(x, y, z, 0) = S 0 (x, y, z), (x, y, z) ∈ G, t = 0, i ∈ M.

Continuous Model Linearization
To obtain sufficient conditions for the existence of a unique solution to the problem (1)-(4), we make additional assumptions about the periodicity of the process: where T > 0-period. Introduce on the surface ∑ of the domain G functions Linearize the source functions in the interval 0 < t < T on a uniform time grid ω τ = {t n = nτ, n = 0, 1, . . . , N; Nτ = T}.
Functions q n i (x, y, z, t n−1 ) are defined at each step of the time grid ω τ . If n = 1, then q 1 i (x, y, z, t 0 ) and it suffices take the functions of the initial conditions q 1 i (x, y, z, 0) ≡ q 0i (x, y, z). When n = 2, . . . , N functions q n i (x, y, z, t n−1 ) = q n−1 i (x, y, z, t n−1 ) are assumed to be known since problem (1)-(4) on the previous time interval t n−2 < t ≤ t n−1 is assumed to be solved.
Then, on the time interval t n−1 < t < t n , Equation (1) has the form: where The initial conditions are attached: where (x, y, z) ∈ G, n = 2, . . . , N, i ∈ M. Border conditions: The source functions divide into two terms: where the term p i (q j ) · q i is linear relative to q i and R q i -nonlinear. Then, the coefficients linear with respect to the q i terms in the functions of the right-hand sides will have the form: , p n−1 Si = 0, and the nonlinear terms will have the form Introduce the operators C and D, which act as follows: Then, the original system can be rewritten as: Introduce an operator L, which acts as follows: Then, system (1) takes the form A Hilbert space H is introduced with the scalar product of vectors x = ξ n 1 , ξ n 2 , . . . , ξ n m , y = η n 1 , η n 2 , . . . , η n m , where n is the time step number, n = 1, N, acting according to the formula: where m is equal to the number of equations of system (1); that is, m = 10.
It is easy to see that this functional for all x ∈ H that satisfies for all the axioms of the scalar product and, therefore, takes place at (Lx, x) > 0. According to [25,26], this condition means the existence of an inverse operator L −1 , which ensures the existence of a solution to problem (1)- (4), and in the case of its positive definiteness, the continuous dependence of the solution on the functions of the right-hand sides and the uniqueness of the solution of the linearized initial-boundary problem.
Find conditions for the positivity of the operator L. To do this, the dot product (Lq, q) > 0 was found, where q = {q i }, i ∈ M, and the corresponding quadratic functional was obtained: According to relations (6), the identities may be written as: Get the functional: Applying the Gauss theorem, the equalities were obtained: Transform relation (11) using formula (12); we arrive at the equality: In accordance with Green's formula and taking into account the boundary conditions (9)-(10), obtain the chain of equalities: Taking into account equality (14), transform the expression for the quadratic functional (13): Let H x , H y , H z be the maximum dimensions of the region G in the horizontal and vertical directions, respectively. The Poincaré inequalities are valid: The expressions on the left-hand side of the functional I are replaced in accordance with the above inequalities by terms not exceeding them, thus constructing the functional I, I ≥ I.
Collect the terms containing q 2 i , i ∈ M: Require that the terms at q n i 2 , i ∈ M be positive at each time level: Write down the obtained conditions for each concentration The coefficients K F j D , K F j E , K PD , K PN , K DN , K 42 , K 23 are positive due to their biological meaning, therefore, only conditions (16), (19)-(21) will be essential. Let us write them down in more detail: The verification of conditions for the existence and uniqueness of a solution to system (1)-(4) is performed layer by layer based on solutions of the initial-boundary value problem.
Using obtained estimations, we are coming to the next theorem.

Theorem 1.
Let an initial-boundary value problem be formulated for a system of equations linearized along the right-hand sides in an interval 0 < t < T on a uniform time grid ω τ = {t n = nτ, n = 0, 1, . . . , N; Nτ = T}, and a chain of initial-boundary value problems (6)-(10) is obtained.
We will assume that as a result of solving these problems, q n−1 i = maxq n−1 i , i ∈ M, are determined based on the solution of the initial-boundary value problem at the previous time level as the maximum solution on the interval t n−1 < t ≤ t n , n = 1, . . . , N.
Let q i belong to the class C 2 (G) the boundary surface of the domain Σ is piecewise smooth, and for each n = 1, N inequality (16), (19)-(21) are satisfied.
Then, the solution to the formulated problem exists and is unique for 0 < t ≤ t 1 and t n−1 ≤ t ≤ t n , n = 2, N.

Numerical Solution of a Multispecies Problem of Phytoplankton Dynamics
To approximate the convective terms in the diffusion-convection-reaction equations, improved Upwind Leapfrog schemes are used, which have better accuracy and a large margin of stability in comparison with those known for large values of the grid Péclet number. To discretize a continuous mathematical model of the dynamics of the most common in Azov Sea phytoplankton summer species (1)-(4), a linear combination of a central difference scheme and an Upwind Leapfrog scheme with weight parameters selected based on minimizing the approximation error were used [29]. The conservativeness of the proposed difference scheme at the discrete level is studied. The conservation of masses (amount of matter) in solving the transfer problem is shown. The fluid volume of the control areas method [30] was used to reduce approximation error at the boundary. The constructed system of discrete equations is implemented using the method of splitting along spatial coordinates (along horizontal and vertical directions). To solve the resulting grid equations, an adaptive modified alternating-triangular iterative method of the variational type is applied, which is a two-layer iterative method with the high convergence rate [31].

Numerical Experiment
Numerical modeling of the problem solution of the phytoplankton populations' dynamics was conducted considering the transformation of the forms of phosphorus, nitrogen, and silicon, using the Azov Sea as an example. A software module was developed based on a mathematical model of biogeochemical cycles, which allows for obtaining threedimensional concentrations distributions of the main phytoplankton populations (green, blue-green, and diatoms) and nutrients (phosphorus, nitrogen, and silicon compounds). The developed module was built into the existing "Azov3D" software package, which allows modeling hydrodynamic processes in the Azov Sea under the influence of winds, the presence of zones with reduced microturbulent exchange in the vertical direction, considering surge phenomena, the Coriolis force, river flows, complex geometry of the computational domain, as well as the rejection of the hydrostatic approximation.
The modeling area corresponds to the physical dimensions of the Azov Sea (355 × 233 km). The size of the grid cell covering the computational area in the horizontal plane is 500 m. The satellite image in Figure 1 visualizes the habitats of green and blue-green algae in the area of the Taganrog Bay and diatoms in the central part of the sea, which are most of the biomass in the warm season according to long-term observations [32][33][34][35][36]. The habitats of phytoplankton populations may change due to changes in the hydrological regime of the reservoir; modeling of such situations is presented in [37].
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 longterm average data for the month of July. Initial concentration values of bluegreen algae (Aphanizomenon flos-aquae)-2.6 mg/L, green algae (Chlorella vulgaris)-2.5 mg/L, diatoms (Sceletonema costatum)-0.9 mg/L; distributions are uniform, with optimal salinity for the first two species phytoplankton-6‰, for the third-12‰, the coefficients of the width of the salinity tolerance interval b 1,2,3 = 2.
The modeling area corresponds to the physical dimensions of the Azov Sea (355 × 233 km). The size of the grid cell covering the computational area in the horizontal plane is 500 m. The satellite image in Figure 1 visualizes the habitats of green and blue-green algae in the area of the Taganrog Bay and diatoms in the central part of the sea, which are most of the biomass in the warm season according to long-term observations [32][33][34][35][36]. The habitats of phytoplankton populations may change due to changes in the hydrological regime of the reservoir; modeling of such situations is presented in [37]. 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 the month of July. Initial concentration values of bluegreen algae (Aphanizomenon flos-aquae)-2.6 mg/L, green algae (Chlorella vulgaris)-2.5 mg/L, diatoms (Sceletonema costatum)-0.9 mg/L; distributions are uniform, with optimal salinity for the  Figure 2 shows the surface distributions of the simulated substances' concentrations: green algae, bluegreen algae, diatoms, phosphates, nitrates, as the most preferred nutrient compounds for all phytoplankton species, and silicon compounds, which are actively consumed by diatoms. first two species phytoplankton-6‰, for the third-12‰, the coefficients of the width of the salinity tolerance interval 1,2,3 2 b = . Figure 2 shows the surface distributions of the simulated substances' concentrations: green algae, bluegreen algae, diatoms, phosphates, nitrates, as the most preferred nutrient compounds for all phytoplankton species, and silicon compounds, which are actively consumed by diatoms. The obtained distributions of substance concentrations are consistent with the data of long-term observations [26,38]. The results of modeling the dynamics of the main phytoplankton populations qualitatively coincide with the data of the Earth's space sensing. The obtained distributions of substance concentrations are consistent with the data of long-term observations [26,38]. The results of modeling the dynamics of the main phytoplankton populations qualitatively coincide with the data of the Earth's space sensing.

Software Implementation
For mathematical modeling of the phytoplankton development dynamics, taking into account the transformation of nutrient forms, a software (SW) was developed. The SW simulates the dynamics of the development of three main types of summer phytoplanktonbluegreen algae (Aphanizomenon flos-aquae), green algae (Chlorella vulgaris), diatoms (Sceletonema costatum); their competition for nutrients, transformation of the forms of these nutrients-phosphorus, nitrogen, and silicon, their consuming, excretion, transition from one biochemical compound to another, the influence of salinity and temperature on the growth rate of phytoplankton are taken into account. The developed software allows modeling biogeochemical processes that determine the biological productivity of such a coastal system as the Azov Sea and the state of the aquatic ecosystem in common [26,37].
In the software module, the calculation took into account the influence of winds, the microturbulent exchange in the vertical direction, considering surge phenomena, the Coriolis force, river flows, complex geometry of the computational domain, as well as the rejection of the hydrostatic approximation. The scheme of the algorithm of the program module "Calculation of the aquatic environment movement" is shown in the Figure 3. SW simulates the dynamics of the development of three main types of summer phytoplankton-bluegreen algae (Aphanizomenon flos-aquae), green algae (Chlorella vulgaris), diatoms (Sceletonema costatum); their competition for nutrients, transformation of the forms of these nutrients-phosphorus, nitrogen, and silicon, their consuming, excretion, transition from one biochemical compound to another, the influence of salinity and temperature on the growth rate of phytoplankton are taken into account. The developed software allows modeling biogeochemical processes that determine the biological productivity of such a coastal system as the Azov Sea and the state of the aquatic ecosystem in common [26,37].
In the software module, the calculation took into account the influence of winds, the microturbulent exchange in the vertical direction, considering surge phenomena, the Coriolis force, river flows, complex geometry of the computational domain, as well as the rejection of the hydrostatic approximation. The scheme of the algorithm of the program module "Calculation of the aquatic environment movement" is shown in the Figure 3. The program block "Calculation of the phytoplankton concentrations" uses a threedimensional water flow velocity vector and takes into account the impact on the phytoplankton development of such abiotic factors as salinity and temperature, the three-dimensional fields of which are obtained as a result of the program module "Calculation of the aquatic environment movement" ("Azov3D.exe").

Results and Discussion
In the field of aquatic ecology Alekseev A.G., Svirezhev Yu.M., Barenblatt I.B., Vinogradov I.M., Medvinsky A.B., Suzuki H., Fukuoka S., Ghosh D., Sarkar P., Jonson B., Ebenman B., Buffoni B., Ruan Yn., and others studied the qualitative properties of analytical models. Their research is based on the classical methods of the theory of differential equations, corresponding to the so-called zero-dimensional models (which do not consider spatial inhomogeneity).
It should be noted that, despite the great attention of scientists to the issues of mathematical ecology, some important problems that serve as a justification for the applicability of a particular model to specific water bodies remain in the shadows. In the field of hydrophysics and biological kinetics, the authors did not find publications concerning the results of an analytical study of the existence and uniqueness of the solution of initialboundary value problems, which are spatially inhomogeneous non-stationary models of biogeochemical cycles for real coastal systems. These issues are important because obtaining the conditions for the existence and uniqueness of solutions, including restrictions on The program block "Calculation of the phytoplankton concentrations" uses a threedimensional water flow velocity vector and takes into account the impact on the phytoplankton development of such abiotic factors as salinity and temperature, the three-dimensional fields of which are obtained as a result of the program module "Calculation of the aquatic environment movement" ("Azov3D.exe").

Results and Discussion
In the field of aquatic ecology Alekseev A.G., Svirezhev Yu.M., Barenblatt I.B., Vinogradov I.M., Medvinsky A.B., Suzuki H., Fukuoka S., Ghosh D., Sarkar P., Jonson B., Ebenman B., Buffoni B., Ruan Yn., and others studied the qualitative properties of analytical models. Their research is based on the classical methods of the theory of differential equations, corresponding to the so-called zero-dimensional models (which do not consider spatial inhomogeneity).
It should be noted that, despite the great attention of scientists to the issues of mathematical ecology, some important problems that serve as a justification for the applicability of a particular model to specific water bodies remain in the shadows. In the field of hy-drophysics and biological kinetics, the authors did not find publications concerning the results of an analytical study of the existence and uniqueness of the solution of initialboundary value problems, which are spatially inhomogeneous non-stationary models of biogeochemical cycles for real coastal systems. These issues are important because obtaining the conditions for the existence and uniqueness of solutions, including restrictions on the coefficients of equations and other input data of the problem, provides specialists with information about the features of the model applicability. As a rule, conclusions on these issues were made not on the basis of theoretical studies, but on the basis of an analysis of approximate solutions during computational experiments for some model problems. The authors of this article tried to fill this gap using the example of a three-dimensional initial-boundary value problem corresponding to a 3D model of geochemical cycles and phytoplankton populations dynamics as applied to the Azov Sea.
To discretize the continuous model, a difference scheme developed by the team of authors was used, which is a linear combination of a central difference scheme and an Upwind Leapfrog scheme. The value of the difference scheme developed by the authors and used in this research is the fact that, in addition to a high order of approximation, this scheme showed its effectiveness in the case of Péclet numbers from 2 to 20, in contrast to the difference scheme with central differences, which is advisable to use for grid Péclet numbers less than 2. For the proposed difference scheme, an increase in the order of the approximation error is allowed, but, in this case, the stencil of the difference scheme increases and ceases to be compact. The central difference scheme for Péclet numbers greater than 2 is unstable.
In the course of previously performed numerical experiments and using data from field measurements of those types of summer phytoplankton concentrations that are included in the model, as well as on the basis of expeditionary studies of the water area, forecasts were obtained for the growing season (April-October) of changes in the number of its individual species, with an error not exceeding 10-15%, acceptable for predicting the processes of hydrophysics and biological kinetics.
The authors of the article performed computational experiments with well-established software systems for solving oceanological problems POM (Princeton Ocean Models), Mars3D, and others with significantly different depths, as well as at critical wind stresses accompanying storm surges. At the same time, the use of models and methods for their numerical implementation and the software Asov3D developed by the authors' team previously allow for the successful reconstruction of the storm surge in the Taganrog Bay on 23-24 September 2014, and earlier, the emergence of a vast zone of hypoxia-anaerobic pollution over an area of more than 1000 square kilometers in the Azov Sea in 2001, as well as the reconstruction of a number of other hazardous phenomena, including eutrophication and "blooming waters", spills of oil and oil products [22][23][24]32].

Conclusions
On the example of the Azov Sea, numerical experiments were carried out on diagnostic and predictive modeling of the processes under consideration on the basis of the developed software package. The modeling results are consistent with the available observational data for the previous period of salinity increase in the Azov Sea in general and the Taganrog Bay in particular. Further research, beyond the scope of this article, will address the emerging trend of significant salinity fluctuations within a 3-4 year observation period and the impact of fluctuations on the species' composition and number of plankton populations.