Currently, special monitoring methods are being developed to observe promptly the changes in the trophic status of water objects [1]. Eutrophication, as a process of uncontrolled increase of primary production in water objects, sometimes can be limited by the nutrients removal from the aquatic ecosystem. In this context, studies aimed at reducing mineral forms of nitrogen and phosphorus consumed by algae are relevant [2]. Mathematical models and simulation are becoming increasingly used tools in water use, for example, in the optimization of wastewater treatment plants [3]. In some cases, advanced mathematical models successfully predict the state of components of aquatic ecosystems [4]. The use of model tools is particularly effective when the calculated values lie within the past range of observations [5].

A sound aquatic ecosystem is in homeostasis [6]. In this case, under the influence of internal regulators, the ecosystem copes with stress and returns to its original state. The adaptability, an important property of the ecosystem, allows it to cope with the consequences of external impact [7]. However, in some cases, in a period of a few years an essential shift of homeostasis indicators occurs (for example, [8]). This is reflected not only in the change in hydrobionts biomasses, but also in the reduction of biodiversity, in replacement of dominant plankton species, etc. [9]. The process of eutrophication is an example of structural changes in the ecosystem of the reservoir [10].

The interaction of phyto- and zooplankton is an important mechanism of the aquatic ecosystem response to climate and anthropogenic changes [11]. Quite a lot of models [12] have been developed in which this relationship is reproduced [13], including specially constructed ones [14]. Still, this mechanism has not yet been properly studied [15]. Most studies report on significant synchronization with a certain time lag of variability of phytoplankton and zooplankton biomass (for example, [16]). In a purpose reseach [17], a positive correlation between these characteristics of hydrobionts was obtained for the George Bank Bay ecosystem. Still, in other cases, such as the Gulf of Maine, there was no stable correlation dependence in biomasses variability for phytoplankton and zooplankton. Friedland et al. [17] note that ‘the inability to form a consistent representation... stresses the need for further research'. It seems that the mechanisms of interaction between phyto- and zooplankton of a real water body cannot be understood without links with other components of the ecosystem.

Simulation of changes in aquatic ecosystem state is necessary to assess its response to the stress factors. In fact, traditional modeling tools (for example, [18], [19]) cannot reproduce long-term trends in ecosystem development required for studying degradation processes. Classical thermodynamic concepts are needed to expand its scope of applicability [20]. A step forward in this direction is the approach of non-equilibrium thermodynamics, which allows to present a structural change as the passage through the chain of conditionally stable states of the aquatic ecosystem. Structural-Dynamic Modeling (SDM) [21] serves as a basis for the study of qualitative changes in the state and assessment of future development stages of aquatic ecosystems.

In this work, we modified the traditional model of the aquatic ecosystem and further used it for estimating the stability loss of the reservoir ecosystem. The trend of long-term changes in the Novosibirsk reservoir eutrophication, the largest reservoir in Western Siberia, was characterized. The interannual variability of components of the reservoir ecosystem was estimated on the basis of modeling for a series of low-water years.

As the basis of the study, the ‘Biogen' model [19], [22] was used, which describes the internal biotic cycle of transformation of nitrogen and phosphorus compounds-the main limiting elements in reservoir ecosystem.

All considered components of the water ecosystem participate in the aforesaid transformation. The mineral forms of nitrogen and phosphorus are consumed by algae and serve as the most important biochemical basis of organically bound substances, which form the biomass of living cells. Along the food trophic chain, organic phosphorus of living matter is included in the cells of organisms of a higher trophic level – zooplankton. As a result of metabolic reactions, dissolved mineral and organic components are released into the surrounding aquatic environment. When the cells of microorganisms die off, a suspended substance, i.e., detritus, which also contains the phosphorus and nitrogen compounds is formed. The role of detritus is important in regeneration of compounds of biogenic elements. During biochemical oxidation, nitrogen and phosphorus bound in detritus are transformed into dissolved organic matter. Thus, the transformation of biogenic compounds in the aquatic environment happens as the complex interaction of aquatic microorganisms and chemicals. Environmental factors influence on concentrations of chemical and biological components in the aquatic ecosystem.

The model of transformation of biogenic compounds, designed to study environmental processes in the reservoir, describes the mechanisms of biochemical transformation of nutrients. Such a model can reproduce the response of the ecosystem under the actual biogenic load and predict probable changes in hydrobiocenosis caused by anthropogenic load decrease or increase in the future. Exchange processes at the ‘bottom sediment-water' boundary are important in the characterized cycle. Therefore, the model describes the transformation of biogenic compounds not only in the aquatic environment, but also in the bottom sediments as well.

The ‘Biogen' model is based on the following fairly strong simplifications. The stoichiometric relationship between the concentrations of carbon (C), nitrogen (N) and phosphorus (P) in the ecosystem components is constant and equal to 106:16:1, respectively. The content of C in the ecosystem accounts for about half of the total abiotic organic mass. The main source of suspended matter in the river is represented by soil particles washed out from the catchment area. Organic matter content in the inflow is 10% of the total concentration of the suspension that corresponds with the organic concentration in the surface layer of soils, prevailing in the watershed. In case of a reservoir, where the main source of suspension is dead plankton, this assumption is unacceptable. Changeability of hydrobionts communities is taken into account through the analysis of dynamics of two aggregated variables: phytoplankton and zooplankton [19], [22].

To study the long-term eutrophication processes of the Novosibirsk reservoir, the following ecosystem components were selected: ammonium nitrogen (N-NH₄), nitrite (N-NO₂), nitrate (N-NO₃), mineral phosphorus (I), dissolved organic matter (C), suspended matter (D), oxygen (O₂). The communities of hydrobionts are represented in the model by two aggregated components: zooplankton biomass (ZO) and phytoplankton biomass (F). The following designations were introduced to describe the interaction of the aquatic environment with bottom sediments: active and passive organic matter (C_B and C_N, respectively), nitrogen and phosphorus compounds sorbed on solid phase (N_S and P_S), interstitial (N_B and P_B). The interaction between ecosystem components is presented in Figure 1.

Figure 1.

A scheme of biochemical components transformation in the aquatic ecosystem according to the model ‘Biogen’

The model equations for estimating average per reservoir volume concentrations (C_i), 1 ≤ i ≤ 9 appear as:

$\frac{d(C_i\times W)}{dt}=W\times R_i+Q^p\times C_i{}^p-Q\_\times C_i+J_i\times \Omega +G_i\times L$ (1)

where i = ZO, F, NH₄, NO₂, NO₃, D, C, I, O₂, W is the reservoir volume, t is the time, R_i is the rate of biochemical transformation of compound C_i, Q^p and $C_i^p$ are the water input discharge and i-th component concentration in the Ob river, Q_– is the water outlet discharge from reservoir, J_i is the mass flow on the interface, Ω is the surface water area of the reservoir, G_i is the lateral loading of diffuse pollution, L is the length of the reservoir shoreline. A list of coefficients with links to calculation methods are presented in the Appendix.

The model ‘Biogen' was calibrated against the observed data collected at the Novosibirsk reservoir, constructed on the Ob river in 1957 (the reservoir is located upper of Novosibirsk city). The reservoir is the largest in Western Siberia. It is used mainly for energy production, navigation and water supply. Its average depth is 8.2 m, the reservoir volume and water-surface areas at the normal affluent level are 8.8 km³ and 1,070 km², respectively, volume of the reservoir drawdown is 4.4 km³. The lake-like part of the reservoir is a site, where the water quality is formed. At normal affluent level, the lake-like part of the reservoir contains 73% of the total reservoir water. With the reservoir drawdown, the share of the lake-like area increases to 90%.

Trophic status of the Novosibirsk reservoir (in the end of 20^th century) was attributed to the oligo-mesotrophic type [23]. The total phytoplankton biomass in the reservoir discharge was insignificant (approximately 0.68 mg/l). In 1981/1982, its maximum (1.1-1.77 g/m³) was observed in late August and September. The species composition of dominating phytoplankton communities underwent numerous changes a year.

A comprehensive study of the Novosibirsk reservoir conducted in 1981/1982 yielded a body of data allowing for a tentative model description of the processes in the reservoir ecosystem [23]. The scenario of annual dynamics of component concentrations in the head section of the reservoir (the Kamen'-na-Obi) was selected based on actual concentrations of mineral forms of N (NH₄, NO₂, NO₃), P, O₂, chlorophyl ‘a', zooplankton, phytoplankton, organic matter and detritus with the use of empirical relationships [24], [25].

Water inflow from the reservoir shore accounted for less than 6% of river runoff [23]. Calculations of N and P fluxes from the reservoir shores, which amount from 0 to 12% of the biogenic load on the river, did not reveal any significant difference for annual variations in the model components. The biogenic load on the reservoir due to inflow from its shores was assumed to be uniformly distributed over months.

Unknown parameters of the ‘Biogen' model were found within the intervals of typical values by minimizing Cr – the Theil's statistical criterion [26] for components of the aquatic ecosystem:

$Cr=\frac{\sqrt{{\displaystyle \sum _{i=1}^n{(X_i-Y_i)}^2}}}{\sqrt{\sum _{i=1}^n{X}_i{}^2}+\sqrt{\sum _{i=1}^n{Y}_i{}^2}}$ (2)

where n is the number of observations, X_i and Y_i are the simulated and observed values of substance concentrations, respectively.

In the 1970s, anthropogenic and climatic changes resulted in eutrophication the Novosibirsk reservoir.

At present, hydrobiologists characterize the reservoir as "a water body of eutrophic type. In accordance with long-term data, average concentrations during the open water period characterize it as weakly eutrophic ..." [27]. Eutrophication often induces almost complete replacement of some species (e.g., phyto- and zooplankton). Rates of these processes differ for various hydrobionts. Simulation of long-term eutrophication, associated with the loss of stable seasonal variability in the ecosystem, requires expanding the boundaries of previous approaches.

The question arises of adequate tools for simulation of long-term succession. It implies the account of changes in species composition and development level of phytoplankton and zooplankton.

The concept of transition to a new trophic status of the ecosystem through a chain of conditionally stable states was used. In this paper, the modification of the ‘Biogen' model proposed in Tskhai and Ageikov [19] is presented. Several variants of sets of key parameters of the model, as opposed to a single one, were considered at each calculation step.

Let us see how the mentioned sets of model parameters were obtained in the work. Tskhai and Ageikov [22] analyzed the sensitivity of the model parameters in advance. The $A_1{}^k$ parameters (Appendix), responsible for the excretory activity of zooplankton (k = 1) and phytoplankton (k = 2), affect the dynamics of aggregated model variables most strongly. Note that a change in the composition of zooplankton and phytoplankton species is a structural change of the model. Functional parameters $A_1{}^k$ (k = 1; 2) of the model vary correspondingly. At each calculation step, three values are selected for k^th parameter:

$\text{(99\%}\times A_1{}^k;A_1{}^k;101\text{\%}{A}_1{}^k)$ (3)

Two aggregated variables of hydrobionts (zooplankton and phytoplankton) are considered in the paper. Therefore, at each calculation step, according to eq. (3), nine variants of the trajectory of the ecosystem model dynamics were simulated. Next, the simulation is performed nine times. The variant with the maximum value of exergy – estimated by eq. (3) – gives the values of $A_1{{}^k}_{\max }$, for the next calculation step. Exergy characterizes the measure of remoteness of the state of a living system from the state of a ‘dead' substance. Thus, an option is selected that increases the adaptability of the living system to environmental conditions. Such a procedure was repeated for the whole simulated period.

A structural change in the ecosystem is a change in the species composition of hydrobionts caused by varying external impacts. The dynamics of a long-term variability of aggregated model variables as a rule depends on structural changes in the ecosystem. The SDM scheme is based on using a set of trajectory variants at each simulated step instead of a single trajectory. To select the optimal trajectory, an optimization problem is solved with the exergy value as the target function. This function application is similar to the Darwin's principle of natural selection as an increase in evolutionary adaptability of ecosystems.

Thus SDM is a simulation on the basis of the heuristic assumptions. An efficiency of the method (as the actual progress in the explanation of ecosystem processes) determines the right of existence of the SDM. The use of the SDM had been possible in more than 23 studies at the time of publication of the review [21]. The results of the next section were obtained through the use of the SDM in the study of long-term trends in the ecosystem state for the Novosibirsk reservoir.

The synergy of techniques for modeling biogeochemical cycles with the non-equilibrium thermodynamics method, described in the previous section, allowed to modify the ‘Biogen' model used in Tskhai and Ageikov [19]. The results of comparison of observed and simulated data on the Novosibirsk reservoir for the ‘Biogen' model (in modified and unmodified forms) are shown in Figure 2 and Figure 3.

Figure 2.

Chemical compounds according to modified and unmodified ‘Biogen' models in the calibration period (1981/1982 hydrological year)

Figure 3.

Biomass of phyto- and zooplankton components according to the modified and unmodified ‘Biogen’ models in the calibration period (1981/1982 hydrological year)

Dashed line ‘A' indicates monthly average observation data. Curve ‘B' presents the result of traditional modeling of biogeochemical cycles. Curve ‘C' is based on the use of the modified model. Numbers on charts are Theil's statistical criteria characterizing the adequacy of simulated and observed results. On each graph, next to ‘B' the value of Cr criterion in traditional modeling is given and next to ‘C' – in modified modeling. The Cr values for both cases satisfy the ecological models [21].

Dynamics in calculated curves obtained from the original and modified models are similar, except for phytoplankton in the summer-autumn period. By and large we deal with a good correlation of calculated and observation data in the both cases.

The Theil's statistical criterion of phytoplankton dynamics decreased, and consequently, the conformity of simulated results and field data improved. Dynamics of O₂, N and P compounds remain practically unchanged when using the modified model unlike unmodified modeling for 1981/1982 year.

The aggregate variable ‘zooplankton' was included to the variables of the model ‘Biogen' [19], [28]. Simulation results of inter-annual variability of zooplankton content are analyzed for the first time in this work.

Actual observations of changes in zooplankton content in 1981 [23] are evidence of two peaks of growth: spring (May-June) and summer-autumn (August) with a period of decrease in total biomass between them (in Figure 3 – the average monthly observation data are marked by stair-step dashed line). This dynamics of biomass is due to division of periods of intensive growth and subsequent death of dominant species of spring and summer-autumn complexes of zooplankton species. In July, the spring zooplankton complex passed the peak of development and the maximum for the summer-autumn zooplankton complex was not apparent for this moment.

When modeling other lakes and reservoirs [12], two peaks in intra-annual variability of total biomass of zooplankton are usually explained by similar trend in phytoplankton biomass. In our case, according to both monthly-averaged observation data and simulation results, the growth peak of the aggregated variable characterizing phytoplankton is single. Spring outburst of diatoms in low water 1981/82 hydrological year was rather short and not intensive. Its contribution did not significantly affect the average monthly biomass of phytoplankton. The difference in variability of phyto- and zooplankton was simulated during low-water year for the ecosystem of the Novosibirsk reservoir for the first time.

In the summer period of low-water year, intense water ‘bloom' occurred due to the dominant algae Aphanizomenon flos-aquae. The scenario of extreme conditions of low-water 1981/82 was repeated during seven years for assessment of the worst-case scenario in terms of water quality. Calculations were made using the modified and non-modified models. In Figure 4, the results of modified model calculation are given.

Figure 4.

Phyto- and zooplankton for seven year simulation with recurring low-water conditions of the calibration period

Of course, the scenario of an exact repetition of conditions for seven years is an idealization because of scarce actual observations. It should be noted that hydrological conditions in the mentioned years were more often low-water ones. Actually the process of eutrophication increase was slower and took not less than 23 years [27]. The result described below is that for the first time, modelling of the sequence of several years for the Novosibirsk reservoir demonstrates a noticeable increase in bloom intensity without increase in external load.

The increase in eutrophication (transition from oligo-mesotrophic to weakly eutrophic), caused by the internal factor (i.e., change in structure) is reproduced through the continuous calculation of several years. This is a novelty, as the unmodified model calculation for the seven identical years gives only a coincidence for whole period.

This suggests that without taking into account structural changes in the ecosystem model, there is no reason to expect the development of eutrophication in the reservoir as a result of the simulation.

The sensitivity analysis of the developed model revealed that excretory activity of hydrobionts plays the key role in biomass change. This factor allows improving the correlation between the calculated and observed data. The mentioned feature of hydrobiont communities is associated with changes in dominant species. It is predicted in Vasiliev [27] that further eutrophication due to peculiar hydrometeorological regime of the reservoir in summer can hardly increase the duration time of cyanoprokaryotic algae bloom, but can surely increase its intensity.

Modified model-based simulation, starting from the second calculated year, shows a continuous increase in the intensity of phyto- and zooplankton development during ‘blooming' in case of a series of low-water years. These calculations not only prove the revealed tendency to increase of intensity in ‘blooming' peak, but also allow us to predict its probable earlier shift.

Predictive calculations for a series of low-water years suggest an annual increase in summer-autumn maximum of zooplankton content at almost constant level of its spring peak. This prediction gives additional arguments in favour of to significant reduction of surplus spills during the period of summer-autumn low water [28].

The accomplished investigation makes it possible to carry out a preliminary assessment of changes in reservoir eutrophication depending on different scenarios of hydrological conditions.

Structural and dynamic modeling of ecosystem processes in this study is based on the concept of exergy. The long-term eutrophication of the Novosibirsk reservoir was simulated for the first time on this basis. For this moment, it is no information about use of this approach based on the concept of exergy for other Russian water bodies. The following natural effects are characterised in this study for the Novosibirsk reservoir for the first time. In accordance with observations in the Novosibirsk reservoir during the 1981, the appearance of two peaks of annual variability of average monthly biomass for zooplankton with single – for phytoplankton was simulated in this study. Also the results obtained suggest a bloom shift in the peak to earlier period of the Novosibirsk reservoir in case of a series of low-water years.