Research of the Solutions Proximity of Linearized and Nonlinear Problems of the Biogeochemical Process Dynamics in Coastal Systems

Abstract

The article considers a non-stationary three-dimensional spatial mathematical model of biological kinetics and geochemical processes with nonlinear coefficients and source functions. Often, the step of analytical study in models of this kind is skipped. The purpose of this work is to fill this gap, which will allow for the application of numerical modeling methods to a model of biogeochemical cycles and a computational experiment that adequately reflects reality. For this model, an initial-boundary value problem is posed and its linearization is carried out; for all the desired functions, their final spatial distributions for the previous time step are used. As a result, a chain of initial-boundary value problems is obtained, connected by initial–final data at each step of the time grid. To obtain inequalities that guarantee the convergence of solutions of a chain of linearized problems to the solution of the original nonlinear problems, the energy method, Gauss’s theorem, Green’s formula, and Poincaré’s inequality are used. The scientific novelty of this work lies in the proof of the convergence of solutions of a chain of linearized problems to the solution of the original nonlinear problems in the norm of the Hilbert space L₂ as the time step τ tends to zero at the rate O(τ).

1. Introduction

Every year, the ecological state of the coastal systems of the South of Russia is deteriorating due to the impact of natural and anthropogenic factors, such as global climate change, increasing water salinity, rising average annual temperatures, agricultural and industrial human activities, etc. Eutrophication of waters, particularly in the Taganrog Bay and the Azov Sea, causes degradation of both individual components of the ecosystem and entire communities of organisms. The collection and analysis of information on the current and predicted states of coastal systems becomes significant. Their variability over time, on a scale of up to several weeks, requires prompt forecasting of adverse phenomena, for which non-stationary spatially inhomogeneous interconnected mathematical models and effective numerical methods for their implementation have already been developed and continue to be improved, allowing us to “play” various scenarios of dynamic biological and geochemical processes in coastal systems.

Coastal systems’ state variability has become the object of many studies by Russian and foreign scientists in the field of biochemistry and biological kinetics. Diagnosis of the state and forecast of changes in the thermo-hydrodynamics of water bodies’ ecosystems is considered in a number of works [1,2]. The article [3] presents estimates of the influence of changes in the hydrological regime on the species composition of phytoplankton and the habitats of the main species. The works [4,5] present the principles of modeling and monitoring the spatio-temporal dynamics of water bodies. A significant amount of work belongs to scientists from the Marchuk Institute of Numerical Mathematics of the Russian Academy of Sciences (INM RAS) and the Shirshov Institute of Oceanology of Russian Academy of Sciences (IO RAS). Research is being conducted on the development of perspective procedures and algorithms for analyzing observational data (ensemble optimal interpolation, Kalman filters, three-dimensional and four-dimensional variational analysis) [6,7]. In works [8,9,10] and others, biogeochemical transformations, hydrochemistry and rates of microbial processes in the water column, the assimilation of chemical elements by hydrobionts of various levels of trophic chains, and the influence of agriculture and industry on the state of the Black, Azov, and Caspian seas are studied. The hydrological regime, models of movement of the aquatic environment, biochemistry, biological kinetics, dynamics of primary bioproduction, and biogenic pollution of water bodies in the South of Russia were studied in [11,12]. The Marine Hydrophysical Institute of Russian Academy of Sciences (MGI RAS) is engaged in the creation and development of a model for diagnosing and forecasting the evolution of the main hydrophysical fields of the Black Sea, which operates in an operational mode (up to 5 days) [13]. Russian and German scientists jointly developed a 1D model of biogeochemical processes [14,15]. This hydrophysical–biogeochemical model simulates the distribution of major nutrients (oxygen, nitrogen, sulfur, phosphorus, manganese, iron) in the pelagic redox layer for the Black and Baltic seas. An example of a three-dimensional model for modeling the processes of hydrodynamics, biological kinetics, and geochemical cycles in coastal systems is MARS [16]. This model is integrated with models of sediment dynamics, microbiology and pollution distribution, and biogeochemical cycles. Cycles of nitrogen transformations in the presence of oxygen and in an anaerobic environment are considered in detail in [17,18]. Extensive studies of the influence of abiotic and biotic factors on phytoplankton communities are described in a number of works by Danish scientists. In [19], the influence of a temperature increase over the course of a month on the development of populations of cyanobacteria and green algae in shallow lakes was studied. In this case, scenarios of low and high nutrient content, nitrogen in particular, were used. The work [20] describes a 1D mathematical model for modeling hydrodynamic and biogeochemical processes in a shallow lake. The modeling of nitrogen cycles in marine systems, including the absorption of nitrogen by bottom sediments, is considered in the work of American scientists [21]. At the same time, various approaches are used to describe denitrification processes: from mechanistic diagenetic models to empirical parameterizations of nitrogen fluxes across the sediment–water interface. A mathematical model of the hydrodynamics of coastal systems, which allows for the modeling of turbulent flows, is described in [22]. The results of the study of the oscillatory behavior of solutions of some classes of differential equations, which allow the modeling of many processes in the ocean, are presented in [23]. Numerical methods and difference schemes used to implement models of hydrodynamics and biogeochemistry are presented in [24,25,26].

This paper presents a three-dimensional non-stationary model of the dynamics of plankton populations and geochemical cycles, which describes the change in the concentrations of the main biogenic substances (compounds of phosphorus, nitrogen and silicon, etc.), considering the mechanisms of plankton development, which allows for the modeling of geographical dynamics of plankton populations, the change in species composition due to changes in biogenic factors, and the development of adverse and dangerous phenomena (rapid flowering of poisonous species of phytoplankton and toxic algae, deadly phenomena). The initial-boundary value problem corresponding to the indicated model contains nonlinear terms, which significantly increases the computational complexity and the time required to solve it. To solve this type of problem, linearization methods are used such as, for example, Newton’s method or, applied to this problem, the method of introducing a time grid with values of nonlinear terms determined with a “delay”, and then solving a system of grid equations using well-established working mechanisms. An overview of modern linearization methods as well as their advantages and disadvantages are presented in [27].

The authors performed a linearization of the initial-boundary value problem, i.e., a replacement of an approximate problem similar to it in terms of dynamic properties, the solution of which is quite close to the solution of the original problem. Linearization is carried out on a time grid; nonlinear coefficients are taken with a delay of one step of the time grid. Next, a problems chain, interrelated by initial conditions and the final solutions, linearized on a uniform time grid is built and, thus, the linearization of the 3D nonlinear model as a whole is carried out.

The use of linearization allows us to replace the functions included in the problem with linear ones and obtain a problem with a linear operator, which can later be solved using numerical methods.

Earlier, in the article [28], the authors’ team presented a description of some aspects of the construction and numerical implementation of the model of biogeochemical cycles and biological kinetics of the multispecies population model. To this end, effective numerical methods and difference schemes have been developed that allow for the simulation of biogeochemical and hydrodynamic processes in a real computational domain of a complex shape [24]. However, some issues of the analytical study of the model remained in the shadows. The authors attempted to close this gap.

It is well-known that the main assumption that allows us to go from a nonlinear problem to a linearized one is the assumption that the deviations of the values of the functions included in the consideration from their analogs, taken as initial ones during linearization, are small. Therefore, questions of “proximity” of solutions to the original nonlinear initial-boundary and linearized problems, on the basis of which the discrete model (difference scheme) was built, are of paramount importance. The scientific novelty of this work lies in the study of the convergence of a chain of linearized problems to the solution of the original nonlinear problems in the norm of the Hilbert space as the parameter, the step of the time grid on which the linearization was carried out, tends to zero. The originality of the work is determined not only by the results obtained, but also by the apparatus used, which is based on the methods of the theory of differential equations.

It is worth noting that the conditions for the existence and uniqueness of the solution of linearized problems were obtained earlier, along with the definition of requirements for the smoothness classes of the input data of the problem [28].

2. Materials and Methods

To model biogeochemical processes in coastal systems, a non-stationary three-dimensional mathematical model based on works [14,18,19], adapted for the Azov Sea, was used. When modeling, the influence of salinity and temperature (abiotic factors) on the growth and death of phytoplankton cells was taken into account. When describing the biogenic-dependent growth rate of phytoplankton, the mathematical equivalent of the Liebig principle was applied using the saturation concentration. The Michaelis–Menten dependence was used to describe the intensity of enzymatic reactions. The model was written on the assumption that in the process of excretion and death, phytoplankton release phosphorus in dissolved and particulate organic forms, which, in the process of phosphatification, turn into an inorganic form—phosphates, which are consumed by phytoplankton. The nitrogen cycle is also described: In the course of life, phytoplankton release nitrogen in an organic form, which decomposes to ammonia. Ammonia in the process of nitrification is oxidized to nitrites, and then to nitrates. The consumption and excretion of silicon by diatoms has also been noted.

2.1. Mathematical Model of Biological Kinetics and Geochemical Cycles

The mathematical model of biological kinetics and geochemical cycles is based on a system of diffusion–convection–reaction with nonlinear coefficients and source functions $R{q}_i$:

$$\frac{\partial q_i}{\partial t}+u\frac{\partial q_i}{\partial x}+v\frac{\partial q_i}{\partial y}+w\frac{\partial q_i}{\partial z}=div\left(\overrightarrow{k}\mathrm{grad}{q}_i\right)+R_{q_i}$$

(1)

where $q_i$ is the concentration of the i-th component, [mg/L]; i∈M, M = {F₁, F₂, F₃, PO₄, POP, DOP, NO₃, NO₂, NH₄, Si}; $\overrightarrow{V}=\left\{u,v,w\right\}$ is the water flow velocity vector, [m/s]; $\overrightarrow{k}=\left(k_h,k_h,k_{\nu }\right)$ is the turbulent exchange coefficient, [m²/s]; and $R_{q_i}$ represents the chemical–biological sources, [mg/(L∙s)]. The system of Equation (1) is written for the following substances: $F_1$ represents the green algae (Chlorella vulgaris), which forms the basis of nutrition for heterotrophs; $F_2$ represents the blue-green algae (Aphanizomenon flos–aquae, cyanobacteria), which competes with green algae for nutrients (toxic); $F_3$ represents the diatom algae (Sceletonema costatum), which consumes silicon, from which it builds frustules—hard shells of cells; $P{O}_4$ represents the phosphates, the mineral form of phosphorus available for consumption by major phytoplankton species; $DOP$ represents the dissolved organic phosphorus, which is released by phytoplankton populations in the process of excretion. In the process of phosphatification, it turns into $P{O}_4$; $POP$ represents the particulate organic phosphorus, which is released by phytoplankton populations in the process of dying off. In the process of autolysis, it turns into $DOP$ and, in the process of phosphatization, into $P{O}_4$; $N{O}_3$ represents the nitrates; $N{O}_2$ represents the nitrites oxidized to nitrates; $N{H}_4$ represents the ammonia oxidized to nitrite. Phytoplankton populations consume all three forms of nitrogen; and $Si$ represents the dissolved inorganic silicon (silicic acids).

Nonlinear functions of the right-hand sides, describing chemical–biological interactions, have the form:

$$R_{F_i}=C_{F_i}(1-K_{F_{i}R})q_{F_i}-K_{F_{i}D}{q}_{F_i}-K_{F_{i}E}{q}_{F_i}, i=\overline{1,3}$$

(2)

$$\begin{gathered} R_{POP}={\displaystyle \sum _{i=1}^3{s}_P{K}_{F_{i}D}{q}_{F_i}}-K_{PD}{q}_{POP}-K_{PN}{q}_{POP} \\ {R}_{DOP}={\displaystyle \sum _{i=1}^3{s}_P{K}_{F_{i}E}{q}_{F_i}}+K_{PD}{q}_{POP}-K_{DN}{q}_{DOP}, R_{P{O}_4}={\displaystyle \sum _{i=1}^3{s}_P{C}_{F_i}\left(K_{F_{i}R}-1\right)q_{F_i}}+K_{PN}{q}_{POP}+K_{DN}{q}_{DOP} \\ {R}_{N{H}_4}={\displaystyle \sum _{i=1}^3{s}_N{C}_{F_i}\left(K_{F_{i}R}-1\right)\frac{f_N^{\left(2\right)}\left(q_{N{H}_4}\right)}{f_N\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)}{q}_{F_i}}+{\displaystyle \sum _{i=1}^3{s}_N\left(K_{F_{i}D}+K_{F_{i}E}\right)q_{F_i}}-K_{42}{q}_{N{H}_4} \\ {R}_{N{O}_2}={\displaystyle \sum _{i=1}^3{s}_N{C}_{F_i}(K_{F_{i}R}-1)\frac{f_N^{\left(1\right)}\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)}{f_N\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)}\cdot \frac{q_{N{O}_2}}{q_{N{O}_2}+q_{N{O}_3}}{q}_{F_i}}+K_{42}{q}_{N{H}_4}-K_{23}{q}_{N{O}_2} \\ {R}_{N{O}_3}={\displaystyle \sum _{i=1}^3{s}_N{C}_{F_i}\left(K_{F_{i}R}-1\right)\frac{f_N^{\left(1\right)}\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)}{f_N\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)}\cdot \frac{q_{N{O}_3}}{q_{N{O}_2}+q_{N{O}_3}}{q}_{F_i}}+K_{23}{q}_{N{O}_2} \\ {R}_{Si}=s_{Si}{C}_{F_3}\left(K_{F_{3}R}-1\right)q_{F_3}+s_{Si}{K}_{F_{3}D}{q}_{F_3} \end{gathered}$$

where $K_{F_{i}R}$ is the specific phytoplankton respiration rate; $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 rate 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; and $s_P$, $s_N$, and $s_{Si}$ are the normalization coefficients between the content of N, P, and Si in organic matter according to the stoichiometric ratio.

The growth rate of phytoplankton is limited by temperature, salinity, and nutrients and is determined by the expressions:

$$\begin{gathered} C_{F_{1,2}}=K_{N{F}_{1,2}}{f}_T\left(T\right)f_S\left(S\right)min\left\{f_P\left(q_{P{O}_4}\right),f_N\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right)\right\} \\ {C}_{F_3}=K_{N{F}_3}{f}_T\left(T\right)f_S\left(S\right)min\left\{f_P\left(q_{P{O}_4}\right),f_N\left(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4}\right),f_{Si}\left(q_{Si}\right)\right\} \end{gathered}$$

where $K_{NF}$ is the maximum specific phytoplankton growth rate.

Functional dependencies on abiogenic factors are described as follows:

$$\begin{gathered} f_T\left(T\right)=\exp \left(-a_l{\left\{\left(T-T_{opt}\right)/T_{opt}\right\}}^2\right), l=\overline{1,3}; f_S\left(S\right)=\exp \left(-b_l{\left\{\left(S-S_{opt}\right)/S_{opt}\right\}}^2\right), l=2,3 \\ {f}_S\left(S\right)=\{\begin{array}{l}{k}_s,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for} \hspace{0.17em}\hspace{0.17em}S\le S_{opt},\\ \exp \left(-b_1{\left\{\left(S-S_{opt}\right)/S_{opt}\right\}}^2\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for} \hspace{0.17em}\hspace{0.17em}S>S_{opt},\end{array} \end{gathered}$$

where $k_s=1$; $T_{opt}$ and $S_{opt}$ are the optimal temperature and salinity for a given type of aquatic organisms; and $a_l>0$ and $b_l>0$ are the coefficients of the width of the range of aquatic organisms’ tolerance to temperature and salinity, respectively.

Functional dependences of the growth rate of phytoplankton on the content of phosphates, silicon, and ammonia are taken in the form of Michaelis–Menten [29], and for nitrates and nitrites, they have the form:

$$f_N^{}(q_{N{O}_3},q_{N{O}_2},q_{N{H}_4})=\left(q_{N{O}_3}+q_{N{O}_2}\right)\exp \left(-K_{psi}{q}_{N{H}_4}\right)/\left(K_{N{O}_3}+\left(q_{N{O}_3}+q_{N{O}_2}\right)\right)$$

where $K_{N{O}_3}$ is the half-saturation constant of nitrates and $K_{psi}$ is the coefficient of ammonium inhibition.

We assume that the coefficients in the functional dependencies describing the sources of the calculated concentrations are positive and independent of time.

Now we consider the corresponding initial and boundary conditions. The area G of the reservoir is a cylindrical area, the side surface of which is formed by the movement of the vertical axis along a piecewise-smooth closed line $\partial {\mathrm{\Sigma }}_H$, and a bounding bottom surface ${\mathrm{\Sigma }}_H$, which is a smooth surface with a “cover” in the chosen coordinate system ${\overline{\mathrm{\Sigma }}}_H={\mathrm{\Sigma }}_H\cup {\mathrm{\Sigma }}_h\hspace{0.17em}$. The bottom base ${\mathrm{\Sigma }}_o$ is the undisturbed free surface of the reservoir, more precisely, the part of it that is “cut out” when the guide (parallel to the Oz axis) moves on the plane Z = 0. We then denote a cylindrical surface enclosed between ${\mathrm{\Sigma }}_H$ and ${\mathrm{\Sigma }}_o$ as $\sigma$. Thus, the closure of the region G, $\overline{G}=G\cup \mathrm{\Gamma }$, where $\mathrm{\Gamma }={\mathrm{\Sigma }}_o\cup {\mathrm{\Sigma }}_H\cup \sigma \hspace{0.17em}$. Considering the introduced notation, the boundary and initial conditions for the system of Equation (1) are formulated as follows:

$$q_i=0\hspace{0.17em}\mathrm{on} \sigma , \mathrm{if} u_{\overrightarrow{n}}<0; \frac{\partial q_i}{\partial n}=0 \mathrm{on} \sigma , \mathrm{if} u_{\overrightarrow{n}}\ge 0$$

(3)

$$\frac{\partial q_i}{\partial z}=0 \mathrm{on} {\mathrm{\Sigma }}_o; \frac{\partial q_i}{\partial z}=-{\epsilon }_i{q}_i \mathrm{on} {\mathrm{\Sigma }}_H$$

(4)

where constants ${\epsilon }_i\ge 0$, i∈M, consider the sinking of phytoplankton to the bottom and their flooding for i∈{F₁, F₂, F₃} and considers the absorption of nutrients by bottom sediments for i∈{PO₄, POP, DOP, NO₃, NO₂, NH₄, Si}.

Given initial values of the components under study, as well as temperature and salinity fields and the water flow velocity vector, at any moment of time for the system of Equation (1):

$$\begin{gathered} q_i\left(x,y,z,0\right)=q_{0i}\left(x,y,z\right), \left(x,y,z\right)\in \overline{G}, t=0 i\in m \\ \overrightarrow{V}\left(x,y,z,0\right)={\overrightarrow{V}}_0\left(x,y,z\right), T\left(x,y,z,0\right)=T_0\left(x,y,z\right), S\left(x,y,z,0\right)=S_0\left(x,y,z\right) \end{gathered}$$

2.2. Continuous Model Linearization

To linearize system (1), we construct a uniform time grid ${\omega }_{\tau }$ with a time step $\tau$ on the time interval $0<t<T$, where T is the characteristic period of biochemical cycles associated with the growing season of the main phytoplankton populations:

$${\omega }_{\tau }=\left\{t_n=n\tau ,n=0,1,\dots ,N;N\tau =T\right\}$$

(5)

At each time step $t_{n-1}<t<t_n$ with numbers $n=1,2,\dots ,N$, we consider the system of Equation (1) linearized with respect to the functions of the right-hand sides $R_{q_i}$, i∈M, the solutions of which are the functions ${\tilde{q}}_{{}_i}^n$, $n=1,2,\dots ,N$, in the form:

$$\frac{\partial {\tilde{q}}_i^n}{\partial t}+\mathrm{div}\left(\overrightarrow{V}\cdot {\tilde{q}}_i^n\right)=\mathrm{div}\left(\overrightarrow{k}\cdot \hspace{0.17em}\mathrm{grad}\hspace{0.17em}{\tilde{q}}_i^n\right)+{\tilde{R}}_{q_i}^n$$

(6)

To Equation (6), we add the initial conditions of the form:

$$\begin{gathered} {\tilde{q}}_i^1\left(x,y,z,0\right)=q_{0i}\left(x,y,z\right), \left(x,y,z\right)\in G \\ {\tilde{q}}_i^n\left(x,y,z,t_{n-1}\right)={\tilde{q}}_i^{n-1}\left(x,y,z,t_{n-1}\right), \left(x,y,z\right)\in G, n=2,\dots ,N \end{gathered}$$

as well as the corresponding boundary conditions, similar to the conditions (3) and (4).

For all t, $t_{n-1}<t<t_n$, $n=1,2,\dots ,N$:

$${\tilde{q}}_i^n=0, \mathrm{if} u_{\overrightarrow{n}}<0; \frac{\partial {\tilde{q}}_i^n}{\partial n}=0, \mathrm{if} u_{\overrightarrow{n}}\ge 0$$

(7)

where $\overrightarrow{n}$ is the outer normal to the side surface $\sigma$, $\left(x,y,z\right)\in \sigma$;

$$\frac{\partial {\tilde{q}}_i^n}{\partial z}=0,\left(x,y,z\right)\in {\mathrm{\Sigma }}_0; \frac{\partial {\tilde{q}}_i^n}{\partial z}=-{\epsilon }_i{\tilde{q}}_i^n, \left(x,y,z\right)\in {\mathrm{\Sigma }}_H$$

(8)

In [28], the functions of the right-hand sides ${\tilde{R}}_{F_i}^n$, $n=1,2,\dots ,N,$ $i\in M$, are defined and a theorem of existence and uniqueness of the solution is formulated under natural constraints on the smoothness of the input data for a linearized problem, with a simplified method of linearization.

3. Investigation of the Proximity of Solutions to the Linearized and Original Initial-Boundary Value Problems by the Energy Method

For the convenience of the study, we formulate the original (nonlinearized) system (1) as a chain of coupled initial-boundary value problems of the form:

$$\frac{\partial }{\partial t}{q}_i^n+\frac{\partial }{\partial x}\left(u{q}_i^n\right)+\frac{\partial }{\partial y}\left(v{q}_i^n\right)+\frac{\partial }{\partial z}\left(w{q}_i^n\right)=\mathrm{div}\left(\overrightarrow{k}\cdot \mathrm{grad}{q}_i^n\right)+R_{q_i^n}^n$$

(9)

where $i\in M$, $\left(x,y,z\right)\in G$, $n=1,2,\dots ,N$, $t_{n-1}<t\le t_n$, and $t\in {\omega }_{\tau }=\left\{t_n=n\tau ,n=1,2,\dots ,N\right\}$, with initial conditions:

$$q_i^1\left(x,y,z,0\right)=q_{0i}\left(x,y,z\right), q_i^n\left(x,y,z,t_{n-1}\right)=q_i^{n-1}\left(x,y,z,t_{n-1}\right), \left(x,y,z\right)\in G$$

as well as boundary conditions (3) and (4) considered on the interval $t_{n-1}<t\le t_n$ for Equation (9).

For brevity, we do not give the definition of functions $R_{q_i}^n$ on each time interval $t_{n-1}<t\le t_n$, considering that:

$$R_{q_i}^n\equiv R_{q_i}\left(x,y,z,t\right), t_{n-1}<t\le t_n$$

Now we consider the specifics of this problem: the specific respiration rate of each of the phytoplankton populations, given by a constant $K_{F_{i}R}$, such that $0<K_{F_{i}R}\ll 1$.

This means that all the functions of the right-hand sides for the biogenic components ($i=4,\dots 10$) $R_i^n$ and ${\tilde{R}}_i^n$ are negative if we focus on the positiveness of all concentrations $q_i^n$ and ${\tilde{q}}_i^n$, $n=1,2,\dots ,N$, and the type of functions $R_i^n$ and ${\tilde{R}}_i^n$, and the typical values of the coefficients, also considering $K_{F_{i}R}-1<0$.

Naturally, such a conclusion is valid if there are no external sources of biogenic concentrations (for example, deficiently treated sewage effluents in a reservoir). We assume that there are initial values of the concentrations of biogenic components $q_i^0\left(x,y,z\right)$ sufficient for the growth of phytoplankton populations in the time interval $0<t\le T^{*}<T$, where $T$ is the considered characteristic period for the entire system (for the Azov Sea, 7–8 months).

Now we turn to the analysis of the closeness of the solutions of the linearized and original problems, considering the assumptions made. We introduce the function:

$$z_i^n\left(x,y,z,t\right)\equiv {\tilde{q}}_i^n\left(x,y,z,t\right)-q_i^n\left(x,y,z,t\right), i\in M, t_{n-1}<t\le t_n, \left(x,y,z\right)\in G,n=1,2,\dots ,N$$

We subtract from Equation (6) each ${\tilde{q}}_i^n$ (number $i$) containing the equation of the number $i$ containing $q_i^n$ from the system of Equation (9), $n=1,2,\dots ,N$. Considering the linearity of the operators participating in the equations, we obtain a system of the form:

$$\frac{\partial z_i^n}{\partial t}+\frac{\partial }{\partial x}\left(u{z}_i^n\right)+\frac{\partial }{\partial y}\left(v{z}_i^n\right)+\frac{\partial }{\partial z}\left(w{z}_i^n\right)=\mathrm{div}\left(\overrightarrow{k}\cdot \mathrm{grad}{z}_i^n\right)+\left({\tilde{R}}_{q_i}^n-R_{q_i}^n\right)$$

where $i\in M$, $n=1,\dots ,N$, $\left(x,y,z\right)\in G$, and $t_{n-1}<t\le t_n$.

It is easy to see that the initial conditions are written as follows:

$$z_i^0=0, i\in M, z_i^n\left(x,y,z,t_{n-1}\right)=z_i^{n-1}\left(x,y,z,t_{n-1}\right), \left(x,y,z\right)\in G, n=1,2,\dots ,N$$

It should be noted that due to the presence of periodicity conditions:

$$z_i^N\left(x,y,z,t_N\right)\equiv 0$$

However, this condition is not used in subsequent calculations.

The boundary conditions are formulated in an obvious way, based on the relation of (7) and (8), as well as that of (3) and (4):

$$z_i^n=0, \mathrm{if} u_{\overrightarrow{n}}<0, \left(x,y,z\right)\in \sigma$$

(10)

$$\frac{\partial z_i^n}{\partial \overrightarrow{n}}=0, \mathrm{if} u_{\overrightarrow{n}}\ge 0, \left(x,y,z\right)\in \sigma$$

(11)

$$\frac{\partial z_i^n}{\partial z}=0, \left(x,y,z\right)\in {\sum }_0$$

(12)

$$\begin{gathered} \frac{\partial z_i^n}{\partial z}=-{\epsilon }_i{z}_i^n, \left(x,y,z\right)\in {\sum }_H \\ i=1,2,\dots ,10, n=1,2,\dots ,N, t_{n-1}<t\le t_n \end{gathered}$$

(13)

We assume that each of the functions $q_i^n$ and ${\tilde{q}}_i^n$ is “square”-integrable in the domain $G$ for all $i\in M$, $n=1,2,\dots ,N$. We introduce the scalar product of functions $\xi \left(x,y,z,t\right)$ and $\eta \left(x,y,z,t\right)$, such that for any $0\le t\le T$, $\left(x,y,z\right)\in G$, there exist bounded integrals ${\displaystyle \underset{G}{\iiint }{\xi }^2\left(x,y,z,t\right)}dxdydz$ and ${\displaystyle \underset{G}{\iiint }{\eta }^2\left(x,y,z,t\right)}dxdydz$, each of which is a continuously differentiable function of the variable $t$.

The scalar product $\left(\xi ,\eta \right)$ is the expression:

$$\left(\xi ,\eta \right)={\displaystyle \underset{G}{\iiint }\xi \left(x,y,z,t\right)}\cdot \eta \left(x,y,z,t\right)dxdydz$$

which is a function depending on the variable $t$. Naturally, we assume that $\mathrm{mes}G=const<+\infty$.

Next, we introduce the Hilbert space $H$ for functions $\xi \left(x,y,z,t\right)$, $\eta \left(x,y,z,t\right)$, …, square-integrable on $G$. We introduce the norm:

$${\Vert \xi \Vert }_{L_2(x,y,z)}\equiv {\left(\xi ,\xi \right)}^{\frac{1}{2}}\equiv {\left({\displaystyle \underset{G}{\iiint }{\xi }^2\left(x,y,z,t\right)dxdydz}\right)}^{\frac{1}{2}}$$

Obviously, each such norm is a non-negative function of the variable $t$, continuously differentiable with respect to this variable. Further, for brevity, the product of differentials $dxdydz$ is denoted as $dG$.

Next, we multiply scalarly both parts of each i-th equation by the function $z_i^n$ and integrate over the variable $t$ from $t_{n-1}$ to $t_n$. We get:

$$\begin{gathered} {\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \frac{\partial z_i^n}{\partial t}dG}\right)dt}+{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \mathrm{div}\left(\overrightarrow{V}\cdot z_i^n\right)dG}\right)dt} \\ -{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \mathrm{div}\left(\overrightarrow{k}\cdot \mathrm{grad}{z}_i^n\right)dG}\right)dt=}{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\left({\tilde{R}}_{q_i}^n-R_{q_i}^n\right)z_i^{n}dG}\right)dt} \end{gathered}$$

(14)

We first transform the left side of the resulting equality, which we denote as $L\left(z_i^n\right)$. We have:

$$\begin{gathered} L\left(z_i^n\right)={\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\frac{1}{2}\cdot \frac{\partial {\left(z_i^n\right)}^2}{\partial t}dG}\right)dt}+{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \mathrm{div}\left(\overrightarrow{V}\cdot z_i^n\right)dG}\right)dt} \\ -{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \mathrm{div}\left(\overrightarrow{k}\cdot grad{z}_i^n\right)dG}\right)dt} \end{gathered}$$

(15)

It is obvious that:

$$\begin{gathered} {\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\frac{1}{2}\cdot \frac{\partial {\left(z_i^n\right)}^2}{\partial t}dG}\right)dt}=\frac{1}{2}{\displaystyle \underset{G}{\iiint }{\left(z_i^n\left(x,y,z,t_n\right)\right)}^{2}dG} \\ -\frac{1}{2}{\displaystyle \underset{G}{\iiint }{\left(z_i^n\left(x,y,z,t_{n-1}\right)\right)}^{2}dG}=\frac{1}{2}{\Vert z_i^n\left(x,y,z,t_n\right)\Vert }_{L_2\left(G\right)}-\frac{1}{2}{\Vert z_i^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{L_2\left(G\right)} \end{gathered}$$

(16)

Here, we take advantage of the fact that:

$$z_i^{n-1}\left(x,y,z,t_{n-1}\right)\equiv z_i^n\left(x,y,z,t_{n-1}\right)$$

Using the Gauss formula (theorem), we transform the second term on the left side of equality (15) considering the boundary conditions (10)–(13). For convenience, we introduce the notation:

$${\mathrm{\Sigma }}^{+}=\{\begin{array}{ll}\left(x,y,z\right)\in \mathrm{\Gamma },& \mathrm{if} \hspace{0.17em}\hspace{0.17em}{u}_{\overrightarrow{n}}^{+}\equiv v_{\overrightarrow{n}}\ge 0\\ \mathrm{\Phi },& \mathrm{if} \hspace{0.17em}\hspace{0.17em}{v}_{\overrightarrow{n}}<0\end{array}$$

We have:

$${\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\cdot \mathrm{div}\left(\overrightarrow{V}\cdot z_i^n\right)dG}\right)dt}={\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }\frac{u_{\overrightarrow{n}}^{+}}{2}}{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }\right)dt}$$

(17)

Next, we transform the expression ${\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{z}_i^n\mathrm{div}\left(\overrightarrow{k}\mathrm{grad}{z}_i^n\right)dG}\right)dt}$, included in the expression for $L\left(z_i^n\right)$. We derive some generalization of Green’s first formula as applied to this expression. Then, we consider a vector function:

$$\overrightarrow{A}\equiv P\cdot \overrightarrow{i}+Q\cdot \overrightarrow{j}+R\cdot \overrightarrow{k}$$

where $P\left(x,y,z,t\right)\equiv z_i^n{k}_h\frac{\partial z_i^n}{\partial x}$, $Q\left(x,y,z,t\right)\equiv z_i^n{k}_n\frac{\partial z_i^n}{\partial y}$, $R\left(x,y,z,t\right)\equiv z_i^n{k}_v\frac{\partial z_i^n}{\partial z}$,$\left(x,y,z\right)\in \overline{G}$, and $t_{n-1}\le t\le t_n$, $n=1,\dots ,N$.

We consider the expression ${\displaystyle \underset{G}{\iiint }\mathrm{div}\overrightarrow{A}\hspace{0.17em}}dG$, which, in accordance with the Ostrogradsky–Gauss formula, can be written as the flow of a vector function $\overrightarrow{A}$ through the boundary surface $\mathrm{\Gamma }$ of the domain $G$:

$${\displaystyle \iiint \mathrm{div}\overrightarrow{A}\hspace{0.17em}dG}={\displaystyle \iint A_{\overrightarrow{n}}}d\mathrm{\Sigma }$$

(18)

where $A_{\overrightarrow{n}}$ is the normal component of vector $\overrightarrow{A}$ (in the direction of the outer normal to $\mathrm{\Gamma }$) on the border $\mathrm{\Gamma }$.

Considering the boundary conditions (10)–(13), it is easy to see that:

$${\displaystyle \underset{Г}{\iint }{A}_{\overrightarrow{n}}d\mathrm{\Sigma }}={\displaystyle \underset{{\mathrm{\Sigma }}_H}{\iint }{k}_v\cdot z_i^n\cdot \frac{\partial z_i^n}{\partial z}}d\mathrm{\Sigma }=-{\displaystyle \underset{{\mathrm{\Sigma }}_H}{\iint }{\epsilon }_i\cdot k_v{\left(z_i^n\right)}^2}d\mathrm{\Sigma }$$

(19)

since the flow is zero on the remaining parts of the boundary $\mathrm{\Gamma }$ (relations (10), (11), (12)).

We then return to the expression ${\displaystyle \underset{G}{\iiint }\mathrm{div}\overrightarrow{A}\hspace{0.17em}}dG$, which can be represented as:

$$\begin{gathered} {\displaystyle \underset{G}{\iiint }\mathrm{div}\left(\overrightarrow{A}\right)}dG={\displaystyle \underset{G}{\iiint }{k}_h{z}_i^n\cdot \left(\frac{{\partial }^2{z}_i^n}{\partial x^2}+\frac{{\partial }^2{z}_i^n}{\partial y^2}\right)}dG \\ +{\displaystyle \underset{G}{\iiint }{z}_i^n\frac{\partial }{\partial z}\left(k_v\frac{\partial z_i^n}{\partial z}\right)}dG+{\displaystyle \underset{G}{\iiint }{k}_h\left[{\left(\frac{\partial z_i}{\partial x}\right)}^2+{\left(\frac{\partial z_i}{\partial y}\right)}^2\right]}dG+{\displaystyle \underset{G}{\iiint }{k}_v{\left(\frac{\partial z_i^n}{\partial z}\right)}^2}dG \end{gathered}$$

(20)

From relations (18), (19), and (20), we have an equality, which can be considered the first Green’s formula in relation to our problem:

$$\begin{gathered} {\displaystyle \underset{G}{\iiint }\left[k_h{z}_i^n\left(\frac{{\partial }^2{z}_i^n}{\partial x^2}+\frac{{\partial }^2{z}_i^n}{\partial y^2}\right)+z_i\frac{\partial }{\partial z}\left(k_v\frac{\partial z_i^n}{\partial z}\right)\right]}dG= \\ =-{{\displaystyle \underset{{\mathrm{\Sigma }}_H}{\iint }{\epsilon }_i{k}_v\left(z_i^n\right)}}^{2}d\mathrm{\Sigma }-{\displaystyle \underset{G}{\iiint }{k}_h\left[{\left(\frac{\partial z_i^n}{\partial x}\right)}^2+{\left(\frac{\partial z_i^n}{\partial y}\right)}^2\right]}dG-{\displaystyle \underset{G}{\iiint }{k}_v{\left(\frac{\partial z_i^n}{\partial z}\right)}^2}dG \end{gathered}$$

(21)

Substituting the right parts of equalities (16), (17), and (21)—the last one with a minus sign—into relation (15), we obtain:

$$\begin{gathered} L\left(z_i^n\right)=\frac{1}{2}{\Vert z_i^n\left(x,y,z,t_n\right)\Vert }_{L_2\left(G\right)}-\frac{1}{2}{\Vert z_i^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{L_2\left(G\right)}+{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }\frac{u_{\overrightarrow{n}}^{+}}{2}{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt \\ +{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{k}_h\left[{\left(\frac{\partial z_i^n}{\partial x}\right)}^2+{\left(\frac{\partial z_i^n}{\partial y}\right)}^2\right]dG+{\displaystyle \underset{G}{\iiint }{k}_v{\left(\frac{\partial z_i^n}{\partial z}\right)}^{2}dG}}\right)dt} \\ +{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }{\epsilon }_i{k}_v{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt \end{gathered}$$

(22)

Let $H_x$, $H_y$, and $H_z$ be the maximum dimensions of the region $G$ in the directions of the coordinate axes $Ox$, $Oy$, and $Oz$, respectively:

$$\begin{gathered} H_x=\underset{\begin{array}{l}{A}^{\prime }\in G\\ A^{″}\in G\end{array}}{\sup }\rho \left(A^{\prime },A^{″}\right), \mathrm{where} A^{\prime }=\left(x^{\prime },y,z\right)\in G, A^{″}=\left(x^{″},y,z\right)\in G \\ {H}_y=\underset{\begin{array}{l}{B}^{\prime }\in G\\ B^{″}\in G\end{array}}{\sup }\rho \left(B^{\prime },B^{″}\right), \mathrm{where} B^{\prime }=\left(x,y^{\prime },z\right)\in G, B^{″}=\left(x,y^{″},z\right)\in G \\ {H}_z=\underset{\begin{array}{l}{C}^{\prime }\in G\\ C^{″}\in G\end{array}}{\sup }\rho \left(C^{\prime },C^{″}\right), \mathrm{where} C^{\prime }=\left(x,y,z^{\prime }\right)\in G, C^{″}=\left(x,y,z^{″}\right)\in G \end{gathered}$$

where $\rho \left(P,Q\right)$ is the Euclidean distance function in $G$. Then, the Poincaré inequalities take place, which are further used to estimate the functional (22) and, as a result, to estimate ${\Vert z_i^n\left(x,y,z,t_n\right)\Vert }_{L_2\left(G\right)}$:

$${\displaystyle \underset{G}{\iiint }{k}_h\left[{\left(\frac{\partial z_i^n}{\partial x}\right)}^2+{\left(\frac{\partial z_i^n}{\partial y}\right)}^2\right]}dG\ge 4k_h\left(\frac{1}{H_x^2}+\frac{1}{H_y^2}\right){{\displaystyle \underset{G}{\iiint }\left(z_i^n\right)}}^{2}dG$$

(23)

$${\displaystyle \underset{G}{\iiint }{k}_v{\left(\frac{\partial z_i^n}{\partial z}\right)}^2}dG\ge 4k_{v_{\min }}\frac{1}{H_z^2}{{\displaystyle \underset{G}{\iiint }\left(z_i^n\right)}}^{2}dG$$

(24)

where $k_{v_{\min }}=\underset{\overline{G}}{\min }\left\{k_v\left(x,y,z\right)\right\}$.

Considering Equations (14) and (22), as well as inequalities (23) and (24), we arrive at the following inequality:

$$\begin{gathered} \frac{1}{2}{\Vert z_i^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le \frac{1}{2}{\Vert z_i^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{{}_{L_2\left(G\right)}}^2 \\ +{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\left(\left({\tilde{R}}_{q_i}^n-R_{q_i}^n\right)z_i^n\right)dG}\right)}dt-{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }\frac{u_{\overrightarrow{n}}^{+}}{2}{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt \\ -{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\left[4k_h\left(\frac{1}{H_x^2}+\frac{1}{H_y^2}\right)+4k_{v_{\min }}\frac{1}{H_z^2}\right]}{\left(z_i^n\right)}^{2}dG\right)}dt \\ -{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }{\epsilon }_i{k}_v{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt \end{gathered}$$

(25)

In what follows, on the right side of the inequality, we omit the terms of the form:

$$-{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }\frac{u_{\overrightarrow{n}}^{+}}{2}{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt \mathrm{and}-{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{{\mathrm{\Sigma }}^{+}}{\iint }{\epsilon }_i{k}_v{\left(z_i^n\right)}^{2}d\mathrm{\Sigma }}\right)}dt$$

Because of this, the right side of inequality (25) only increases. To complete the evaluation, we consider the expression:

$${\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }\left(\left({\tilde{R}}_{q_i}^n-R_{q_i}^n\right)z_i^n\right)dG}\right)}dt\equiv S_i^n$$

For simplicity, we consider the case $i=F_1$.

Based on equalities (2) and (5), we conclude that the value of this expression is maximum when there is no limit on the growth of the phytoplankton population by nutrients and abiogenic factors, that is, the population growth rate is the highest. Then, the value of the expression $S_1^n$ is the largest and is also the largest (”worst”) estimate of the error. Based on the maximum value of the coefficient $C_{F_1}^n$, which is equal to the constant $K_{N{F}_1}$, we come to the estimate:

$$S_1^n\le \left[K_{N{F}_1}\left(1-K_{F_{1}R}\right)-K_{F_{1}D}-K_{F_{1}E}\right]\cdot {\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{\left(z_i^n\right)}^2}dG\right)}dt$$

(26)

We note that, in accordance with the hydrobiological meaning of the constants included in the expression in square brackets, $0<K_{F_{1}R}\ll 1$, $K_{F_{1}D}>0$, and $K_{F_{1}E}>0$ and up to a certain point in time, when $C_{F_1}^n=K_{N{F}_1}$, the concentration may increase. Considering inequalities (25) and (26), we obtain (for $i=F_1$),

$$\begin{gathered} {\Vert z_1^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le {\Vert z_1^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{{}_{L_2\left(G\right)}}^2+2[K_{N{F}_1}\left(1-K_{F_{1}R}\right)-K_{F_{1}D}-K_{F_{1}E} \\ -4\left[k_h\left(\frac{1}{H_x^2}+\frac{1}{H_y^2}\right)+k_{v_{\min }}\frac{1}{H_z^2}\right]]\cdot {\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{\left(z_1^n\left(x,y,z,t\right)\right)}^2}dG\right)}dt \end{gathered}$$

(27)

We note that:

$$z_1^n\left(x,y,z,t_n\right)=z_1^{n-1}\left(x,y,z,t_{n-1}\right)+\frac{\partial z_1^{n-1}}{\partial t}\left(x,y,z,t_{n-1}^{*}\right)\left(t-t_{n-1}\right)$$

(28)

where $t_{n-1}\le t_{n-1}^{*}\le t_n$, $t_{n-1}<t_{n-1}^{*}\le t$, and $t-t_{n-1}\le \tau$.

Since all partial derivatives $\frac{\partial {\tilde{q}}_i^{n-1}}{\partial t}$ and $\frac{\partial q_i^{n-1}}{\partial t}$ are assumed to be continuous at the closure of the domain of definition, they are bounded; therefore:

$$\left|\frac{\partial z_1^{n-1}}{\partial t}\right|=\left|\frac{\partial {\tilde{q}}_i^{n-1}}{\partial t}-\frac{\partial q_i^{n-1}}{\partial t}\right|\le K_1^{n-1}\equiv const>0$$

(29)

Therefore, from (28) and (29), we have:

$$\begin{gathered} {\left(z_1^n\left(x,y,z,t\right)\right)}^2\le {\left(z_1^{n-1}\left(x,y,z,t_{n-1}\right)+K_1^{n-1}\tau \right)}^2\le \\ \le 2\cdot \left({\left(z_1^{n-1}\right)}^2+{\left(K_1^{n-1}\right)}^2{\tau }^2\right) \end{gathered}$$

(30)

Then we have, considering (30):

$$\begin{gathered} {\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}\left({\displaystyle \underset{G}{\iiint }{\left(z_1^n\left(x,y,z,t\right)\right)}^2}dG\right)}dt={\displaystyle \underset{G}{\iiint }\left({\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}{\left(z_1^n\left(x,y,z,t\right)\right)}^{2}dt}\right)}dG\le \\ \le 2\tau {\displaystyle \underset{G}{\iiint }{\left(z_1^{n-1}\left(x,y,z,t_{n-1}\right)\right)}^2}dG+2{\left(K_1^{n-1}\right)}^2{\tau }^3\mathrm{mes}G \end{gathered}$$

(31)

where $\mathrm{mes}G\equiv {\displaystyle \underset{G}{\iiint }dG}$.

Considering estimates (27) and (31), we have the inequality:

$$\begin{gathered} {\Vert z_1^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le (1+4\tau [K_{N{F}_1}\left(1-K_{F_{1}R}\right)-K_{F_{1}E}-K_{F_{1}D} \\ -4\left[k_h\left(\frac{1}{H_x^2}+\frac{1}{H_y^2}\right)+k_{v_{\min }}\frac{1}{H_z^2}\right]])\cdot {\Vert z_1^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{{}_{L_2\left(G\right)}}^2 \\ +2{\left(K_1^{n-1}\right)}^2{\tau }^3\mathrm{mes}G \end{gathered}$$

(32)

For convenience, we denote:

$$\begin{gathered} \alpha \tau \equiv 4\tau \left[K_{N{F}_1}\left(1-K_{F_{1}R}\right)-K_{F_{1}E}-K_{F_{1}D}-4\left[k_h\left(\frac{1}{H_x^2}+\frac{1}{H_y^2}\right)+k_{v_{\min }}\frac{1}{H_z^2}\right]\right] \\ \beta {\tau }^3\ge 2{\left(K_1^{n-1}\right)}^2{\tau }^3\mathrm{mes}G, \beta =\left(\underset{1\le n\le N}{\max }{\left(K_1^{n-1}\right)}^2\right)\cdot 2\mathrm{mes}G \end{gathered}$$

Then, estimate (32) can be written in a compact form as follows:

$${\Vert z_1^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le \left(1-\alpha \tau \right){\Vert z_1^{n-1}\left(x,y,z,t_{n-1}\right)\Vert }_{{}_{L_2\left(G\right)}}^2+\beta {\tau }^3$$

It is easy to distinguish (can be proved by induction) that the chain of inequalities–equalities is true:

$$\begin{gathered} {\Vert z_1\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le \left(1+\alpha \tau \right)\left({\Vert z_1^{n-2}\left(x,y,z,t_{n-2}\right)\Vert }_{{}_{L_2\left(G\right)}}^2\left(1+\alpha \tau \right)+\beta {\tau }^3\right) \\ +\beta {\tau }^3={\left(1+\alpha \tau \right)}^2{\Vert z_1^{n-2}\left(x,y,z,t_{n-2}\right)\Vert }_{{}_{L_2\left(G\right)}}^2+\left(1+\alpha \tau \right)\beta {\tau }^3+\beta {\tau }^3\le \\ \le {\left(1+\alpha \tau \right)}^2\left({\Vert z_1^{n-3}\left(x,y,z,t_{n-3}\right)\Vert }_{{}_{L_2\left(G\right)}}^2\left(1+\alpha \tau \right)+\beta {\tau }^3\right)+\left(1+\alpha \tau \right)\beta {\tau }^3+\beta {\tau }^3= \\ ={\left(1+\alpha \tau \right)}^3{\Vert z_1^{n-3}\left(x,y,z,t_{n-3}\right)\Vert }_{{}_{L_2\left(G\right)}}^2+\left({\left(1+\alpha \tau \right)}^2+\left(1+\alpha \tau \right)+1\right)\beta {\tau }^3\le \\ \le {\left(1+\alpha \tau \right)}^n{\Vert z_1^0\left(x,y,z,t_0\right)\Vert }_{{}_{L_2\left(G\right)}}^2+\beta {\tau }^3{\displaystyle \sum _{K=1}^n{\left(1+\alpha \tau \right)}^{K-1}} \end{gathered}$$

(33)

We note that:

$$z_1^0\left(x,y,z,t_0\right)\equiv 0$$

(34)

and the sum $S_n\equiv {\displaystyle \sum _{K=1}^n{\left(1+\alpha \tau \right)}^{K-1}}$ is estimated as follows: It is obvious that

$$\begin{gathered} S_n=\frac{{\left(1+\alpha \tau \right)}^{n-1}}{{\left(1+\alpha \tau \right)}^{n-1}}{\displaystyle \sum _{K=1}^n{\left(1+\alpha \tau \right)}^{K-1}}= \\ ={\left(1+\alpha \tau \right)}^{n-1}\left(1+\frac{1}{1+\alpha \tau }+\frac{1}{{\left(1+\alpha \tau \right)}^2}+\dots +\frac{1}{{\left(1+\alpha \tau \right)}^{n-1}}\right) \end{gathered}$$

(35)

The final sum in parentheses on the right side of equality (35) is the sum of a finite number (n terms) of a convergent geometric progression, in which the first term is 1 and the denominator is $q\equiv \frac{1}{1+\alpha \tau }$. Then, the sum is written as:

$$\begin{gathered} S_n={\left(1+\alpha \tau \right)}^{n-1}\left(1+\frac{1}{1+\alpha \tau }+\dots +\frac{1}{{\left(1+\alpha \tau \right)}^{n-1}}\right)= \\ =\left(\left(1-\frac{1}{{\left(1+\alpha \tau \right)}^n}\right){\left(1+\alpha \tau \right)}^{n-1}\right)/1-\frac{1}{1+\alpha \tau }=\frac{{\left(1+\alpha \tau \right)}^n-1}{\alpha \tau }\le \frac{{\left(1+\alpha \tau \right)}^n}{\alpha \tau }\le C_0\frac{e^{\alpha T}}{\alpha \tau } \end{gathered}$$

(36)

where $C_0=const$. Considering (33), (34), and (36), we obtain the following estimate:

$${\Vert z_1^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2\le C_0\beta \frac{{\tau }^3}{\alpha \tau }\le \frac{C_0\beta }{\alpha }{\tau }^2$$

The last inequality implies an estimate that guarantees the closeness (convergence at $\tau \to 0$) of the solutions of the linearized and nonlinear problems for the substation $F_1$ (initial function $q_{F_1}\equiv q_1$) in $L_2\left(G\right)$ on the sequence of grids ${\omega }_{\tau }\left(\tau \to 0\right)$, and the fulfillment in the case of any $\tau$ inequality:

$$\underset{C_1\equiv const>0}{{\Vert z_1^n\left(x,y,z,t_n\right)\Vert }_{{}_{L_2\left(G\right)}}^2}\le \underset{n=1,2,\dots ,N}{C_1\tau }$$

It should be noted that from the condition of periodicity $q_{F_1}\left(x,y,z,0\right)=q_{F_1}\left(x,y,z,T\right)$ and the chosen method of linearization of the problem, it follows that $\Vert z_1^N\left(x,y,z,t_N\right)\Vert =0$.

In addition, there is a (perhaps not the only) moment of time $T_1^{*}$, $T_1^{*}<T$, after which the decrease in the error begins. Non-uniqueness can be caused by the entry of biogenic components into area $G$ from any source. By analogous reasoning, one can prove that:

$${\Vert z_2^n\left(x,y,z,t_n\right)\Vert }_{L_2\left(G\right)}\le C_2\tau$$

where $C_2$ and $C_3$ are positive constants.

4. Discussion

Analytical study of a mathematical model is an important stage of mathematical modeling, and should precede computational experiments. Often, this step is skipped and researchers immediately proceed to the application of numerical methods and computational experiments. Previously, in [29], the authors’ team researched the continuous problem (1) with initial and boundary conditions attached to it. For the proposed model, nonlinear source functions are linearized on a uniform time grid, wherein the values of the nonlinear terms are determined as their final values in the previous time layer (with a delay) and a chain of initial and final solutions of diffusion–convection–reactions initial-boundary value problems is formed. Sufficient conditions for the solutions’ existence and uniqueness to initial-boundary value problems of the dynamics of plankton populations and biogeochemical cycles for a system of partial differential equations in the Hilbert space L₂ are determined. For this, quadratic functionals are constructed and, on the basis of the energy method, the Gauss and Poincaré theorems are referred to, inequalities are determined that guarantee their positivity, and the solutions existence, uniqueness and continuous dependence on the right-hand side functions. In this work, the analytical study of the model is continued. Along with the study of the correctness of the formulation of the initial-boundary problem of biological kinetics and geochemical processes, the convergence of the above chain of linearized problems to the solution of the initial nonlinear problems in the L₂ space norm was carried out as the time step tended to zero. In the course of the study, it was shown that the sequence of solutions of the resulting chain of linearized equations, as the time step tends to zero, tends to solve the original nonlinear problem in the L₂ space norm.

According to the authors, there are no similar results in this area. The analytical study of the described model of biological kinetics and geochemical cycles is quite new, since the authors of other works were focused on practical results and forecasting. There are practically no analytical studies of the correctness and uniqueness of the solution of the formulated problems.

5. Conclusions

The authors studied the mathematical aspects of the model of geochemical cycles and biological kinetics of a multispecies model of populations, considering the following factors: the movement of the water flow, 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 species of plankton populations, including their growth, reproduction, natural decrease in numbers, etc. The input data for the model of biogeochemical processes are of the form of the velocity vector of the aquatic environment, which is calculated on the basis of the model of hydrodynamics of coastal systems [30]. Combining these two models allows us to simulate the researched processes in a medium with a significant density gradient and a large difference in depths and considers the complex shape of the computational domain. The authors’ team developed difference schemes that allow for the consideration of these features of coastal systems [24].

The analytical study presented in this article allowed us to apply numerical methods to the proposed model of biogeochemical processes and to conduct a number of computational experiments. For example, the geographic dynamics of phytoplankton populations in the northeastern part of the Azov Sea (Taganrog Bay) was studied depending on changes in the hydrological regime [27], and forecasts were made for changes in algal habitats under conditions of increased water salinity. The results of diagnostic modeling are in good agreement with field data; the obtained distributions of phytoplankton concentrations qualitatively coincide with satellite images. A further goal of the study of this model is to determine, using numerical analysis, a set of parameters for problems of the dynamics of plankton populations and biogeochemical cycles, on the basis of the obtained analytical estimates, the infinitesimal change of which leads to bifurcation of the solutions of the dynamic system. This would explain such phenomenon as the emergence of stable spatial structures of phytoplankton. We also plan to validate the constructed model based on Earth remote sensing data.