Ecological state evaluation for Lake Glubokoe in Moscow region

This study aims to assess the trophic state of a small lake in modern climatic conditions using a one-dimensional model called MyLake, developed at the Norwegian Institute for Water Research. This temperate Lake Glubokoe, situated in the western part of Moscow region, Russia, is chosen as the study object. It is a small but deep dimictic lake, classified in publications as mesotrophic. Its hydrological and water quality observations were mostly conducted in the first half of the 20th century and are only occasional in later years. The model is calibrated on extensive hydrological observation data collected during the growing season of 2017. It is validated on observation data of 1991–2015. After the model is adapted via an automatic optimization, it is able to accurately describe the lake’s thermal regime and dynamics of the parameters needed for the trophic state evaluation – the mineral phosphorus and chlorophyll concentrations. The Lake Glubokoe trophic state dynamics, based on the average summer chlorophyll concentration in the surface layer, is reviewed from a model simulation of its hydro-ecological regime over 1991–2015. The mean trophic state over this period is determined as β-mesotrophic, reaching an α-eutrophic state in periods of intensive algal growth. These estimates correlate well with the observations of the zooplankton dynamics in the lake.


Introduction
The climate change and anthropogenic impact on aquatic ecosystems are causing significant acceleration of their eutrophication. A warmer climate increases the warming of the water column and intensifies the density stratification during summer [1]. At the same time, the rate of hypolimnion oxygen consumption also increases, while anoxic conditions are increasingly arising, leading to the emergence of biogenic substances from bottom sediments. When thermal stratification is upset, these substances get into the epilimnion, where they become available for phytoplankton consumptionwhich leads to more frequent outbreaks of "blooming".
Lake systems are much more vulnerable to pollution than reservoirs, due to the former's longer water cycle. It is very hard to reverse lake degradation, and lake recovery, if at all possible, may take decades. At the same time, the self-cleaning potential of small lakes is lower than that of medium and large ones, and this particular type of water body is the most common one in Russia (there are more than 2.5 million small lakes in the country).
Studies of long-term changes in the ecological state of aquatic ecosystems revealed increased eutrophication in numerous lakes located in various zones, manifested in increased bacterial activity (which affects the rate of hypolimnion oxygen consumption) and increased phytoplankton biomass 2 and chlorophyll-a concentrations [2][3]. Similar trends were discovered even in lakes located in nature reserves, i.e., unaffected by anthropogenic pollution [4]. Predictive calculations show that if the climate keeps getting warmer, blooming of toxic blue-green algae in lakes will increase even further [5][6].
The fifth evaluation report of the Intergovernmental Panel on Climate Change (IPCC) points out that climate warming was noted in Russia since the 1950s, and at the beginning of the 21st century its impact on the environment only increased [7][8]. At the same time, in recent decades the network of water quality monitoring stations on Russian rivers, lakes, and reservoirs was steadily contracting. Due to insufficiently frequent observations, it is becoming practically impossible to reliably assess the state of aquatic ecosystems based on field data. Meanwhile, the assessments of the lakes' trophic state made several decades ago no longer reflect the current situation. Under such circumstances, hydroecological modeling would seem to be an effective research techniquean advanced tool for studying lake water quality and forecasting further transformation of aquatic ecosystems due to the climate change or anthropogenic impact.
At this stage, it is practically impossible to distinguish between the impact of natural factors causing changes in the trophic state of lakes and that of anthropogenic ones, since the ecosystems of the vast majority of water bodies are transformed by human activity. The impact of the climatic factor specifically can be accurately assessed only for water bodies whose ecosystems have not been affected by anthropogenic transformation. The number of such water bodies is small and continues to rapidly decline, especially in Moscow Region where the level of urbanization has been skyrocketing in recent decades.
Accordingly, we have assessed the current trophic state of one small lake located in the Central European area of Russia. The long-term dynamics of its ecosystem was affected almost solely by natural factors. Due to the lack of recent observation data the assessment was based on mathematical modeling.

Site description
Lake Glubokoe (55.75 N 36.51 E) is located in the European Russia in the western part of Moscow Region within the boundaries of a state natural reserve and, thus, its catchment area was only marginally affected by human activities. The lake is surrounded by mixed spruce-broadleaf forests growing on loam soils, which cover about 80% of its catchment area; the shores are heavily waterlogged.
The N. Zograf hydrobiological station is located on the lake shore; it is one of the first limnological stations in the world and the oldest one in Russia. The staff has been monitoring the lake's zooplankton community since the station's establishment, but the water quality in the lake has not been observed for more than 50 years. Since the second half of the 20th century, the lake's water chemistry dynamics was studied only sporadically.
The shape of the lake basin is close to conical. The maximum length is 1,200 m, the width is 850 m, the maximum depth reaches 32 m, and the average depth is 9.3 m. The surface area is 593,000 m 2 , the volume is about 5.5 million m 3 , and the length of the shoreline is 3,414 m. Water from the catchment area comes mainly from swamps and ditches. The precipitation provides 74% of the total inflow and, thus, the lake's water has low mineral content. The bottom of the lake in its deepest part is covered with silt which mostly consists of detritus, while near the shores the bottom sediments consist of sand, clay, and peat.
The annual water level amplitude is 40-60 cm, and the total variability reaches 1.1 m. In terms of the mixing regime Lake Glubokoe can be classified as dimictic. Freezing begins in October-December and thawing in April-May. Stratification, which happens during summer, leads to a rapid decrease in the hypolimnion oxygen content and the formation of anoxia. The average water temperature in the surface layer during the ice-free period, according to observations for the last 30 years, ranges from  21.2 °C (2010). The average annual water temperature in the surface layer over the entire observation period was about 18.8° C. The lake has low biological productivity; the average transparency is 3.8 m. In terms of the trophic state index (TSI), the lake belongs to the mesotrophic type [9]. The long-term monitoring of the crustacean zooplankton community and the coastal monitoring of the main physical and chemical water quality indicators revealed a tendency towards eutrophication in recent years. However, due to lack of regular observations of the phytoplankton community growth, it is not possible to confirm or disprove this hypothesis without modeling.

Data
In 2017 we performed hydro-ecological monitoring of Lake Glubokoe in order to assess its current state and obtain detailed data, which subsequently were used to calibrate the model. The monitoring activities included measuring the vertical distribution of the water temperature, the electrical conductivity, the dissolved oxygen content and pH, and sampling at different depths to determine the content of mineral and total phosphorus, nitrates, ammonium, silicon, the concentration of major ions, the content of organic substances, and the concentration of chlorophyll-a. Observations were carried out throughout the vegetation period once every two weeks.
Since 1991 the hydrobiological station monitors the water level, the water transparency, the water surface temperature, and the pH level in the coastal zone. This allowed verifying the model using an array of data covering a long-term period.
The meteorological data used in this paper were obtained at the Mozhaysk weather station (55.52 N 36.00 E) located 40 km to the south-east of the object of research. Mean daily values of the air temperature, total cloudiness, relative humidity, wind speed, and the daily rainfall were used as the model inputs.

MyLake model
One-dimensional models are best suited for studying small lakes. In view of the fact that in small lakes vertical heterogeneity in the distribution of characteristics is much more significant than the horizontal 4 one, averaging out of the layers' characteristics seems to be acceptable. Therefore, it becomes possible to schematically view the lake as a single vertical. One-dimensional models are widely used to study the thermal regime and eutrophication of lakes.
Recently, process-based models became more popular than empirical ones. Since the former are based on fundamental equations, they can be used to study a wide range of lakes located in very different landscapes.
The one-dimensional MyLake model we chose as the study tool was designed by the Norwegian Institute for Water Research (NIVA) [10] and successfully applied to study lakes in Norway, Finland, Canada, the USA, etc. [11][12][13][14]. Typically MyLake is employed to measure the impact of climate change on thermal and ice regimes and study algae blooming outbreaks [15][16]. The model's code is written for the MATLAB environment; it is free for use and open for user modifications. We used the model version 1.2.1.
MyLake allows reproducing the average daily vertical distribution of water temperature, the concentration of phosphorus compounds, phytoplankton (from the chlorophyll concentration), and ice cover dynamics. The latest versions of the model allow one to study the dynamics of dissolved oxygen and dissolved organic carbon [17].
Variables such as layer area, depth, and diffusion coefficient are estimated at layer interfaces, while the volume and most of the model variables (temperature, the concentration of phosphorus and chlorophyll-a) are averaged over each layer. The thickness of layers is set by the user and may vary between 0.1 m and 3 m.
The main equations applied in the MyLake model to calculate the temperature profile are the thermal conductivity equation and the equations of potential and kinetic energy balance; the results provide the basis for assessing the intensity and depth of wind mixing. If the kinetic energy of the wind impact is higher than the potential energy in the layer (calculated from the density profile), the epilimnion is mixed with the layer located below it, and the accumulated potential energy is subtracted from the kinetic energy. The procedure is repeated until the remaining kinetic energy becomes lower than the potential one.
The advection-diffusion equation, which essentially is a heat conductivity equation with an added advective term, is used to calculate the dynamics of dissolved and suspended substances.
The phytoplankton dynamics is modeled on the basis of a classical scheme in which the phytoplankton growth rate depends on the concentration of nutrients, the water temperature, and the solar radiation. The phytoplankton death rate is estimated from the reaction of first-order chemical kinetics.
The meteorological data required for the modeling include solar radiation, average daily air temperature, total precipitation, relative humidity, air pressure, and wind speed. The MyLake model does not actually calculate radiation and turbulent heat fluxes; the formulae for their calculation from observation data were included via a special utility MATLAB Air-Sea Toolbox. It implements standard bulk expressions for the calculation of latent and sensible heat fluxes and wind stress as a simplified version of a TOGA/COARE code by Fairall et al [18]. The net long-wave radiative flux is determined using the formula of Berliand [19] with correction for cloud cover. Also, if the input data contain no solar radiation, the flux of short-wave radiation on the lake surface is calculated.
The inputs also include the initial distribution of water temperature, total and mineral phosphorus, and chlorophyll-a content. The lake morphometry is taken into account by setting the bathymetric curve.  [10]. The quality of water temperature simulations largely depends on the vertical diffusion coefficient and the wind sheltering coefficient. The more sensitive ecological parameters include settling velocity for chlorophyll, the half-saturation level for inorganic phosphorus, the growth and mortality rates of phytoplankton, settling velocity for chlorophyll, the PAR saturation level for photosynthesis, etc. (Table 1).

Calibration and validation of MyLake model
The model parameters were optimized automatically using the modFit function of the FME package [20] in the R language, which uses the Nelder-Mead algorithm to find the local extrema (minima) for the function [21].
The upper and lower limits of the variation range, along with the initial values, were set for each parameter based on published data. After that the function identified the vector of parameters which provided the best convergence of the modeled and actual data; the standard root mean square error (RMSE) was chosen as the model quality criterion during optimization. Calibration of the parameters of the physical and ecological modules was carried out simultaneously and, thus, the cumulative relative change in individual variables RMSE was set as the objective function. The optimization procedure was repeated several times under different initial conditions to obtain robust estimates.
After optimization, the convergence of the modeled results with observation data was estimated by calculating the following statistics: the mean absolute error (MAE), the root mean square error (RMSE), and the Theil inequality coefficient (U) [22], which is widely applied in ecological modeling (the lower the U value, the better the convergence between the observed and modeled data): . Optimisation of the model parameters allowed us not only to minimize the RMSE, but also achieve good convergence of the characteristics of vertical distribution. The model adequately reproduces the lake's thermal structure and its transformation during the vegetation period: warming of the epilimnion, formation and deepening of the thermocline, and the autumn cooling of the water column. The vertical temperature distribution during summer is the main factor that affects the chemical and biological processes and, thus, high quality of the thermal regime modeling is necessary for studying the ecological state of lakes. Calibration of the chlorophyll-a content parameters also allowed us to achieve good convergence of the modeled results with the observation data: in most cases the calculated diagrams match the actual ones ( Figure 2). When assessing the quality of the calculation it should also be taken into account that with a low content of photosynthetic pigments the error of the spectrophotometric method of 6 determining chlorophyll-a concentration is 10%, and at concentrations less than 1 μg/L it is 20%. Note that the RMSE for the calculation of the chlorophyll-a concentration is 1.83 µg/L, and the Theil coefficient is 0.024, and so the result can be considered good. The model performed least successfully in reproducing the dynamics of mineral phosphorus compounds. A probable cause of the errors is lack of data on the phosphorus input into the lake from the lateral inflow. Based on the modeling results, the phosphorus content near the surface is close to zero due to its consumption by phytoplankton, whereas according to the observation data this does not occur. However, as in the case of chlorophyll, it should be kept in mind that the error margin of the laboratory procedure for determining the content of phosphates in water is 5%, and the lower limit of sensitivity is 5 μg/L, while the concentrations in the analysed samples were lower than this value, and 7 so they cannot be considered fully reliable. The model reproduces the dynamics of phosphorus distribution satisfactorily: the RMSE value is 3,62 µg/L, and the Theil inequality coefficient is 0,061 ( Table 2).
The obtained RMSE, MAE, and Theil inequality coefficient values indicate high quality of the models of the lake ecological state. When the calibration of the model parameters was completed, verification was carried out for the period from 1991 to 2015 based on observations of the surface water temperature (Fig.3). The results of the verification for 1991-2015 showed a mean square error of 1.85 °C, a mean absolute error of -0.044 °C, and a Theil coefficient value of 0.047. The correlation between the modeled and observed water temperature values has the approximation coefficient R 2 = 0.922, with the length of the series N = 262, which indicates high quality of the simulated thermal regime of the lake. A detailed analysis of the water temperature dynamics during the ice-free period in specific years with different meteorological conditions reveals a match between the dynamics of the actual and calculated temperatures.

Assessment of the lake ecological state
The results of the modeling of the ecological regime of Lake Glubokoe were used to analyze the dynamics of the lake trophic status, which was determined from the average concentration of chlorophyll-a in the surface layer during the summer period. The average content of chlorophyll-a in the surface layer of Lake Glubokoe between June and August for the period from 1991 to 2015 was 11.7 μg/L, which, according to the trophic state index values, allows us to classify the lake as mesotrophic. The maximum values in the surface layer during the summer months range between 16.3 and 23.3 μg/L and, on average, amount to 19.5 μg/L, which corresponds to the eutrophic type (Fig. 4). The highest peak concentrations of chlorophyll do not necessarily occur in years with higher average summer values, and vice versa. For example, in the summer of 1999 the average chlorophyll content was low, but the peak concentration of 22.9 μg/Lthe second highest value over the simulation periodwas recorded. A probable explanation is that the highest daily values can be caused by brief bursts of blooming combined with low photosynthetic activity of phytoplankton during the whole summer. The highest peak chlorophyll-a concentrations in the surface layer of the lake were simulated in 1999, 2000, and 2003, at 23.3, 22.9, and 22.4 μg/L, respectively. The highest mean summer concentrations of photosynthetic pigments during summer were recorded in 1994, 2000, and 2015. The model data on the summer average chlorophyll concentrations for the period in question  were compared with the meteorological data for Moscow region published in [22]. The summers of 2002 and 2010 were extremely hot (the number of days with peak temperatures exceeding 25 °С was 45 and 64, respectively). However, no highest peak and summer average chlorophyll-a content values were observed during these years. On the contrary, the average summer chlorophyll-a concentrations in 2002 and 2010 were the lowest ones (10.7 and 9.6 μg/L, respectively). This was probably due to excessive insolation, and the fact that extremely high water temperature does not promote phytoplankton growth. Atmospheric droughts were also observed in 2010; the amount of precipitation during summer was only 59 mm. The low input of nutrients required for phytoplankton growth (mineral phosphorus) with rainfall resulted in low chlorophyll-a concentration. In the years with moderately hot summers (30-40 days with peak temperatures exceeding 25 °С), the chlorophyll content in the surface layer was higher than in 2002 and 2010. In such years short periods of cooling promote deepening of the thermocline and increased input of phosphorus from the hypolimnion, where the anoxic conditions stimulate its recovery from the bottom sediments. Also, such years are favorable for more active photosynthesis, since the phytoplankton growth is not limited by excessive insolation.
The chlorophyll content in the years with moderate temperature during the vegetation period (13-29 days with peak temperatures over 25 °С) was close to the average long-term value.
In recent decades the trophic status of Lake Glubokoe has not been evaluated due to lack of observations of the chlorophyll concentration and phosphorus content. Therefore, modeling was the only way to perform a detailed study of the year-to-year changes in the lake ecology.

Conclusions
The use of mathematical modeling in the course of this study allowed us to obtain data describing the current and recent hydro-ecological state of Lake Glubokoe. The model adaptation and automatic optimization of the model parameters, combined with verification, helped us achieve a high-quality simulation of the lake thermal dynamics and accurately describe the trophic state of the lake based on the calculated mineral phosphorus and chlorophyll content values. According to the model results, the trophic state of the lake from 1991 to 2015 on average was estimated as mesotrophic, reaching an eutrophic state during periods of maximum algal growth. In general, based on analyzing the long-term dynamics of the chlorophyll concentration during the last 30 years, it can be concluded that there is a trend towards eutrophication of the lake. The model estimates of the trophic status generally coincide with the results of long-term monitoring of zooplankton in the lake [23].