Energy Budget, Water Quality Parameters and Primary Production Modeling in Lake Volvi in Northern Greece

Abstract

A lake’s water quality and its ecosystem structure are mainly determined by heat storage change through its energy budget, dissolved oxygen, nutrients and primary productivity. A one-dimensional water quality model for lakes was used to estimate temperature, dissolved oxygen, phytoplankton (as chlorophyll-α), and inorganic and organic phosphorus. Evaporation, energy budget and surface water temperature were also computed. The results of the mathematical model simulation are presented and evaluated. Data from Lake Volvi in Central Macedonia (in northern Greece) for three successive years (2013 to 2015) were used to calibrate and recalibrate the model. The model results of water temperature, dissolved oxygen and primary productivity (Chlα) were compared with measurements for the years 2013 to 2015. The comparison showed that the predicted values of these parameters were all in good agreement with the measurements. The simulation results of water quality parameters generally exhibited the same seasonal dynamic and inter-annual variations as the measured data. The simulation results of the model application provided important information on changes in the physical, chemical and biological variables of the lake. The water temperature and heat fluxes at the water–atmosphere interface are crucial variables related to climate changes.

1. Introduction

Many biotic and abiotic factors affect and determine lake water quality, with temperature, primary productivity, dissolved oxygen and nutrient concentrations the most significant of them.

The Water Framework Directive of the European Community (2000/60) and the relevant law of Greece suggest the use of mathematical models in water resources management in order to predict quantitative and qualitative changes. The results of mathematical models are used to overcome the weaknesses of in situ measurements and sampling and, in some cases, to find the most suitable sampling site. These models can be important tools in the management of a lake’s aquatic system and useful tools for management to quantify the potential effects of climate change, such as the deterioration in ecological quality through increased Chlα concentration, increased frequency of heatwaves [1], and to estimate the emission of greenhouse gases from inland water bodies [2,3].

Mathematical models have been widely used for studying the water temperature and quality of lakes during the last decades [4]. They can connect hydrodynamics with chemical and biological processes. The eddy diffusion (or differential) models have been successfully used to simulate the distribution of thermal stratification, temperature and water quality parameters in a variety of limnological studies [5,6,7,8].

Bonnet and Poulin (2004) [9] developed the one-dimensional model DyLEM-1D, which can describe the seasonal variation of different species of phytoplankton and the cycles of main nutrients and oxygen in lakes and reservoirs. Arhonditsis and Brett (2005) [10] created a model to simulate eutrophication in lake Washington in the USA and to check alternative management scenarios.

The WASP model [11] was developed by the U.S. EPA and is a mathematical model to analyze water quality problems in water resource systems and can be applied in one, two or three dimensions. The CAEDYM model [12] is a biochemical mathematical model of water quality, which was designed to be used in conjunction with hydrodynamic models, such as the DYRESM and ELCOM models. The DYRESM model [13] is a one-dimensional model of vertical temperature distribution in lakes and reservoirs. In addition, the DYRESM-WQ model [14] can also describe water quality parameters.

The models CE-QUAL-W2 [15] and CE-QUAL-R1 [16] were developed by the U.S. Army Corps of Engineers. CE-QUAL-W2 is a two-dimensional (vertical and transverse) hydrodynamic model which can be applied in lakes, reservoirs, rivers, and rivers’ deltas. CE-QUAL-R1 is a one-dimensional mathematical model describing the vertical distribution of temperature and the chemical and biological parameters in reservoirs. Mi et al. (2020) [17] used the model CE-QUAL-W2 to simulate the metalimnetic oxygen minimum in the Rappbode Reservoir (Germany). Model MINLAKE2012, a lake thermal modeling approach, was used to evaluate the potential climate impacts on the thermal structure and biology of lakes [18]. Xu et al. (2022) [19] combined the LHS and Morris uncertainty and sensitivity analysis methods to examine the uncertainty and sensitivity of main water quality parameters in Tai Lake in China.

The development and use of mathematical models to study the quality, quantity, and biological variables of the lakes in Greece are very sparse and limited. Babajimopoulos and Papadopoulos (1986) [5] developed and applied a one-dimensional model to predict the vertical distribution of water temperature in Lake Vegoritis. Antonopoulos and Gianniou (2003) [20] and Gianniou and Antonopoulos (2007, 2014) [21,22] developed and presented a one-dimensional mathematical model which can simulate temperature, dissolved oxygen, phosphorus and phytoplankton distribution, as well as the water evaporation in lakes. The model was applied successfully in Lake Vegoritis in Greece.

The study of quality, eutrophication, and the management and protection of many lakes, coastal lagoons and artificial reservoirs in Greece using mathematical models of eutrophication and water quality has been a subject of discussion during the last decades [23,24,25,26,27]. The trophic state and the nutrient loads of lake Karla [28] and the water resources management in the Strymon river and Lake Kerkini [29] have been investigated in more recent years.

Artificial neural network (ANN) models are an alternative to the above process-based deterministic models, where multiple physical and chemical factors and meteorological variables are used to estimate thermal, water quality and eutrophication components. They have been widely used in hydrological and environmental processes [30,31]. Different machine learning models have been used to forecast lake surface water temperature [32,33]. Other artificial neural network models have been used for chlorophyll-α prediction [34,35].

One of the most important water impoundments in Greece is Lake Volvi, the second largest lake (in size) in the country. It is located in a hydrological basin with the twin lake Koronia (or Lagadas), whose overflow ran to Lake Volvi until the middle of the 1990s. Significant environmental pressures have been received since the 1990s due to heavy industrial and agricultural pollution and water resource mismanagement in the catchment area [36].

The aim of this work is the calibration and verification of the QUALAKE model, which describes the energy (heat balance and evaporation), water quality (temperature, dissolved oxygen, nutrients—mainly inorganic and organic phosphorus) and biological productivity (phytoplankton as chlorophyll-α) in lakes. The calibration and verification of the QUALAKE model in different lake ecosystems increase the value and credibility of the model. The model QUALAKE was used to simulate these processes in Lake Volvi in Central Macedonia in Greece. The structure of the model was improved in some parts to better describe the different processes, and the results of the water quality variable simulations are presented. QUALAKE can be a useful tool for studying the energy and water quality status of the lake and for estimating the critical external loading to Lake Volvi, at which the lake will meet EU Water Framework Directive (WFD) requirements. The only suitable data for calibrating the model during the last decade was the limited available data collected from the National Water Monitoring Network from 2012 to 2015 [37].

2. Material and Methods

2.1. Site Description

Lake Volvi, along with Lake Koronia, is located in northern Greece and occupies the deeper part of the hydrological basin of Mygdonia (40°37’ N, 23°21’ E) at an altitude of 37 m above sea level (in 2012). Lake Koronia overflowed into Lake Volvi through a ditch until the 1990s. Lake Volvi overflows into the Strymonikos gulf through the Richios stream. The surface elevation, volume, surface area, and maximum depth of the lake in 2014 were 37.0 m above sea level, 1000 × 10⁶ m³, 75.0 km², and 23.5 m, respectively. The trophic status of the lake is considered oligotrophic to mesotrophic.

The lakes have a combined hydrological basin of a total area of 2024 km², with the Lake Volvi subbasin occupying 1247 km². Figure 1 shows the site in Greece, the hydrological basin, with the two subbasins of lakes Volvi and Koronia and the bathymetric map of Lake Volvi.

The lakes Volvi and Koronia are considered and protected by the Ramsar Convention as a system of water bodies of international importance for the value of the wetland habitat, ideal for a variety of flora and fauna species [38].

2.2. Water Surface Energy Budget

The water surface energy budget [32,39,40] is described by:

Q_t = R_n − H − LE

(1)

where Q_t is the change in the energy (thermal) content of the water body (MJ m⁻² day⁻¹), R_n is the net radiation (MJ m⁻² day⁻¹), H is the energy conducted to or from the body of water as sensible heat (MJ m⁻² day⁻¹), and LE is the latent heat flux (due to evaporation) (MJ m⁻² day⁻¹).

The net radiation (R_n) is described as:

$${\mathrm{R}}_{\mathrm{n}}={\mathrm{R}}_{\mathrm{s}}-{\mathrm{R}}_{\mathrm{s}\mathrm{r}}+{\mathrm{R}}_{\mathrm{a}}-{\mathrm{R}}_{\mathrm{a}\mathrm{r}}-{\mathrm{R}}_{\mathrm{b}\mathrm{r}}=(1-{\alpha }_{\mathrm{s}}){\mathrm{R}}_{\mathrm{s}}+(1-{\alpha }_{\mathrm{a}}){\mathsf{\epsilon }}_{\mathrm{a}}\mathsf{\sigma }{({\mathrm{T}}_{\mathrm{a}}+273)}^4-\mathsf{\epsilon }\mathsf{\sigma }{({\mathrm{T}}_{\mathrm{s}}+273)}^4$$

(2)

where R_s and R_a are the short-wave radiation incident to water surface and the incoming long-wave radiation from the atmosphere, respectively (MJ m⁻² day⁻¹); R_sr and R_ar are the reflected short-wave and long-wave radiation, respectively (MJ m⁻² day⁻¹); R_br is the back (long-wave) radiation emitted from the body of water (MJ m⁻² day⁻¹); T_a is the air temperature above the water surface (°C); T_s is the surface water temperature (°C); α_s is the reflectivity of the sort-wave radiation of the water surface (α_s = 0.07); α_a is the reflectivity of the long-wave radiation of the water surface (α_a = 0.03); ε_a and ε are the atmospheric and water emissivity, respectively; and σ is the Stefan–Boltzmann constant (σ = 4.89 × 10⁻⁹ MJ m⁻² day⁻¹ K⁻⁴).

The flux of sensible heat (H) is related to the evaporative heat flux (LE) through the Bowen ratio (β) [41]:

$$\mathsf{\beta }=\frac{\mathrm{H}}{\mathrm{L}\mathrm{E}}=\gamma \frac{{\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{a}}}{{\mathrm{e}}_{\mathrm{s}\mathrm{w}}-{\mathrm{e}}_{\mathrm{d}}}$$

(3)

where e_sw is the saturation vapor pressure at the temperature of the water surface (kPa), e_d is the air vapor pressure above the water surface (kPa), and γ is the psychrometric constant (in kPa °C⁻¹).

The change in the thermal content of the water body (Q_t) is calculated using the change in the lake’s temperature for the time step [21,22], according to the equation:

$$Q_{\mathrm{t}}=\frac{{\rho }_{\mathrm{w}}{c}_{\mathrm{p}\mathrm{w}}}{A_{\mathrm{s}}}\frac{\mathrm{d}}{\mathrm{d}t}{\displaystyle \underset{0}{\overset{{\mathrm{z}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\int }}AT\mathrm{d}\mathrm{z}}$$

(4)

in which c_pw is the specific heat of water (ΜJ kg⁻¹ °C⁻¹), A_s is the surface area of the lake (m²), A(z) is the horizontal area as a function of depth (m²), and T(z,t) is the water temperature (°C) as a function of depth (z) and time (t).

2.3. Mathematical Modeling

The one-dimensional, eddy diffusion model QUALAKE has been used in the past [21,22] to study the water quality, evaporation and energy budget of Lake Vegoritis in Greece. The model uses the method of finite elements to solve the differential equations which describe the spatial and temporal changes of temperature and concentrations of phytoplankton (as the concentration of chlorophyll-α), soluble reactive phosphorus (SRP), organic phosphorus (OP) and dissolved oxygen (DO). The model was calibrated and verified using measured data of water quality parameters for different years [20,21,22]. The model results, during calibration and recalibration, showed that there was good agreement between simulated and measured values of the examined parameters [22].

The description of equations for the simulation of temperature, Chlα, phosphorus and DO distributions in the lake have been presented elsewhere [21,22]. In the following paragraphs, the equations for the submodels for water temperature, SRP, OP, primary production (expressed as the production of phytoplankton biomass) and dissolved oxygen are presented.

2.3.1. Lake Water Temperature Modeling

The heat transport is based on the one-dimensional vertical diffusion equation of the form [7,41,42]:

$$\frac{\partial \mathrm{T}}{\partial \mathrm{t}}=\frac{1}{\mathrm{A}}\frac{\partial }{\partial \mathrm{z}}(\mathrm{K}\mathrm{A}\frac{\partial \mathrm{T}}{\partial \mathrm{z}})-\frac{1}{\mathrm{A}}\frac{\partial (\mathrm{A}\mathrm{q})/\partial \mathrm{z}}{{\mathsf{\rho }}_{\mathrm{w}}{\mathrm{c}}_{\mathrm{p}\mathrm{w}}}$$

(5)

where K is the eddy diffusion coefficient (m² day⁻¹), q(z) is the internal distribution of heat sources due to solar radiation absorption inside the water column (MJ m⁻² day⁻¹), ρ_w is the water density (kg m⁻³), c_pw is the specific heat of water (MJ kg⁻¹ °C⁻¹), A(z) is the horizontal area as a function of depth (m²), and T(z,t) is the water temperature (°C).

The internal distribution of solar radiation q(z) is given by [4]:

$$\mathrm{q}(\mathrm{z})=(1-{\mathsf{\beta }}_{\mathrm{s}}){\mathrm{q}}_{\mathrm{s}\mathrm{n}}\exp (-\mathsf{\eta }\mathrm{z})$$

(6)

where η is the extinction coefficient for solar radiation in water (m⁻¹) and β_s is the fraction of net short-wave radiation absorbed at the surface of the lake.

The eddy diffusion coefficient (K) is calculated using the following equation [20,43]:

$$\mathrm{K}={\mathrm{K}}_0\mathrm{f}$$

(7)

where K₀ is a reference value of K (m² day⁻¹), and f is a function of the Richardson number. The Richardson number depends on the vertical gradient of water temperature, the friction velocity (w*) and the coefficient of expansion. The value of K₀ is related to the friction velocity, which is a function of wind velocity [41,44]. Details of the parameter f as a function of the Richardson number are given in Sundaram and Rehm (1973) [43] and Henderson-Sellers (1984) [41].

2.3.2. Water Quality Parameters Modeling

The one-dimensional vertical diffusion of each of the water quality parameters (SRP, OP, Chlα and DO) is described by a diffusion equation of the following form [9,45,46]:

$$\frac{\partial {\mathrm{C}}_{\mathrm{k}}}{\partial \mathrm{t}}=\frac{1}{\mathrm{A}}\frac{\partial }{\partial \mathrm{z}}(\mathrm{A}{\mathrm{K}}_{\mathrm{D}\mathrm{C}}\frac{\partial {\mathrm{C}}_{\mathrm{k}}}{\partial \mathrm{z}})+{\mathrm{S}}_{\mathrm{k}}$$

(8)

in which C_k is the concentration of the k water quality parameter k (C_k equals C_SRP, C_OP, C_Chlα, C_DO, respectively) (mg m⁻³) as a function of depth (z, in m) and time (t, in days); A is the horizontal area of the lake at the depth z (m²); K_DC is the vertical diffusion coefficient (m² day⁻¹); and S_k is a term describing the net production rate of the k parameter in each layer (mg m⁻³ day⁻¹).

The diffusion coefficient of each water quality parameter is directly related to the vertical thermal stratification structure [46,47,48] and is usually identical to the thermal diffusivity coefficient, which is the sum of the coefficients of molecular diffusivity (D_o, in m² day⁻¹) and turbulent diffusivity (E(z,t), in m² day⁻¹). It is estimated in each time step using a number of equations as presented by Antonopoulos and Gianniou (2003) [20].

The term S_k for the diffusion equation of SRP (k = SRP) includes the losses of SRP due to phytoplankton uptake and water outflows, as well as its sources from phytoplankton respiration, mineralization of organic phosphorus and water inflows in each layer, according to the equation:

$${\mathrm{S}}_{\mathrm{S}\mathrm{R}\mathrm{P}}=-{\mathrm{a}}_{\mathrm{p}\mathrm{c}}\mathrm{P}\mathrm{H}+{\mathrm{a}}_{\mathrm{p}\mathrm{c}}\mathrm{R}+\mathrm{M}\mathrm{I}{\mathrm{N}}_{\mathrm{O}\mathrm{P}}+\frac{1}{\mathrm{A}\mathrm{\Delta }\mathrm{z}}({\mathrm{Q}}_{\mathrm{i}}{\mathrm{C}}_{\mathrm{S}\mathrm{R}\mathrm{P}}{}_{\mathrm{i}}-{\mathrm{Q}}_0{\mathrm{C}}_{\mathrm{S}\mathrm{R}\mathrm{P}})$$

(9)

where a_pc is the phosphorus to chlorophyll-α concentration ratio in phytoplankton cells; PH is the phytoplankton growth rate (mg m⁻³ day⁻¹); R is the phytoplankton respiration loss (mg m⁻³ day⁻¹); Δz is the vertical spatial step (m); MIN_OP is the term describing the mineralization of OP to SRP (mg m⁻³ day⁻¹); Q_i and Q₀ are the inflows and outflows of water from the lake, respectively (m³ day⁻¹); and C_SRPi and C_SRP are the concentrations of SRP in inflow water and in the lake, respectively (mg m⁻³).

The term S_k for the diffusion equation of OP (k = OP) includes the losses of OP due to mineralization and the gains from the mortality of organisms in the lake, the water inflows and outflows and the losses due to sinking of OP:

$${\mathrm{S}}_{\mathrm{O}\mathrm{P}}={\mathrm{a}}_{\mathrm{p}\mathrm{c}}\mathrm{M}-\mathrm{M}\mathrm{I}{\mathrm{N}}_{\mathrm{O}\mathrm{P}}+\frac{1}{\mathrm{A}\mathrm{\Delta }\mathrm{z}}({\mathrm{Q}}_{\mathrm{i}}{\mathrm{C}}_{\mathrm{O}\mathrm{P}}{}_{\mathrm{i}}-{\mathrm{Q}}_0{\mathrm{C}}_{\mathrm{O}\mathrm{P}})-\frac{1}{\mathrm{A}}\frac{\partial (\mathrm{A}{\mathrm{w}}_{{\mathrm{s}}_{\mathrm{O}\mathrm{P}}}{\mathrm{C}}_{\mathrm{O}\mathrm{P}})}{\partial \mathrm{z}}$$

(10)

where M is the mortality rate of phytoplankton (mg Chlα m⁻³ day⁻¹) including grazing by zooplankton, w_sop is the sinking velocity of OP (m day⁻¹), C_OPi is the concentration of OP in the inflowing water (mg m⁻³), and the other terms are already defined.

The term S_k for the diffusion equation of Chlα (k = Chlα) includes phytoplankton production by photosynthesis and the losses of phytoplankton biomass due to respiration, mortality and zooplankton grazing and the sinking of Chlα to sediments:

$${\mathrm{S}}_{\mathrm{C}\mathrm{h}\mathrm{l}\mathrm{a}}=\mathrm{P}\mathrm{H}-\mathrm{R}-\mathrm{M}-\mathrm{G}\mathrm{Z}-\frac{1}{\mathrm{A}}\frac{\partial (\mathrm{A}{\mathrm{w}}_{\mathrm{s}\mathrm{C}\mathrm{h}\mathrm{l}\mathrm{a}}{\mathrm{C}}_{\mathrm{C}\mathrm{h}\mathrm{l}\mathrm{a}})}{\partial \mathrm{z}}$$

(11)

where GZ is the loss due to zooplankton grazing (mg m⁻³ day⁻¹) and w_sChlα is the vertical velocity of phytoplankton sedimentation (sinking velocity) (m day⁻¹).

The term S_k for the diffusion equation of DO (k = DO) includes oxygen production by photosynthesis and the losses due to respiration and the biochemical oxygen demand:

$${\mathrm{S}}_{\mathrm{D}\mathrm{O}}={\mathrm{\alpha }}_{\mathrm{o}\mathrm{c}}\mathrm{P}\mathrm{H}-{\mathrm{\alpha }}_{\mathrm{o}\mathrm{c}}\mathrm{R}-{\mathrm{S}}_{\mathrm{B}\mathrm{O}\mathrm{D}}$$

(12)

where α_oc is the oxygen to chlorophyll-α concentration ratio and S_BOD is the loss of biochemical oxygen demand (mg m⁻³ day⁻¹).

The boundary condition at the water surface (z = 0) for SRP, OP and Chlα diffusion equations is considered to be zero flux across the boundary [22,48], while the boundary condition for DO is considered to be flux type due to the reaeration processes [21,22].

The bottom boundary condition (z = b) for SRP is defined by a flux term due to the release of SRP from the bottom sediments [21,22], while the bottom boundary condition (z = b) for the DO equation is described by the oxygen flux to satisfy the sediment oxygen demand (SOD). The bottom boundary condition (z = b) for OP and Chlα is defined by an equation of zero flux along the boundary [48].

2.4. Results Evaluation

The quantitative evaluation of the model results is based on the statistical criteria of root mean square error (RMSE) and correlation coefficient (r²) between the measured and simulated values [49]. The optimum value of RMSE criteria is zero, while r² is one.

2.5. Data for Simulations

The meteorological variables for the study area were measured at a station located at the north-west of the lake at a distance of 35 km (National Observatory of Athens, Weather station of Lagadas) at 40°42′ N, 23°06′ E with an elevation of 87 m. These measured data include the daily values of maximum, minimum and average air temperature (T_max,T_min, and T_ave, respectively), relative humidity (RH_av), wind velocity (u₂) and precipitation (Pr).

The climate of the study area is considered as semi-arid Mediterranean. The mean annual precipitation ranges between 450 and 600 mm, with a mean annual rain depth of 519 mm. The mean annual temperature is 13.3 °C, with the maximum monthly temperature occurring in July and the minimum in January.

Table 1. Basic metrological data for the three years of simulation and for the time period from 2013 to 2019.

		Tave, °C	RHave, %	Rainfall, mm	u2, m s−1	Rs, MJ m−2day−1
2013	Ave	15.48	75.5	420.60	0.624	15.04
	Max	29.3	97	25.8	2.641	28.37
	Min	−2.5	44	0	0.000	0.76
2014	Ave	15.68	75.5	914.70	0.827	14.08
	Max	29.1	97	118	3.420	29.19
	Min	−2.5	44	0	0.000	0.01
2015	Ave	14.37	71.5	573.60	0.463	15.29
	Max	28.8	98	59.6	4.222	27.2
	Min	−4.1	47	0	0.000	0.31
2013–2019	Ave	15.01	71.77	539.33	0.372	15.54
	Max	30.4	98	1024.7	3.683	30.02
	Min	−8.3	39.5	396.6	0.000	0.01

presents the maximum, minimum, and average values of temperature (°C), relative humidity (%), wind speed (m s⁻¹), solar radiation (MJ m⁻² day⁻¹) and rainfall for each year of simulation and the average values for the time period from 2013 to 2019.

In the absence of sunshine duration or solar radiation data, the multi-linear regression (MLR) equation of solar radiation to meteorological variables from the meteorological station of the Aristotle University of Thessaloniki (AUTH), which is located north-west at a distance of 42 km (40°54′ N; 23°06′ E), was used [50]. This multi-linear regression (MLR) equation uses the meteorological variables T_max, T_min, T_av, RH_av, R_a, TD (=T_max − T_min) and u₂, with the following forms:

$${\mathrm{R}}_{\mathrm{s}}=11.0267+0{.47768\mathrm{R}}_{\mathrm{a}}+0.83564\left(\mathrm{TD}\right)-2{.2386\left(\mathrm{TD}\right)}^{1/2}-0{.16562\mathrm{RH}}_{\mathrm{av}}$$

(13)

The coefficient of determination (r²) of the model was 0.846, based on the daily dataset for the years 2011–2013 from the AUTH station.

The annual water balance of Lake Volvi, estimated with the MIKE HYDRO model using historical data from the period 1970 to 2000 [51], consists of an average direct precipitation of 421 mm, evaporation of 729 mm, runoff of 282 mm, pumping from the lake of 9.33 mm, outflow to the aquifer of 8.8 mm, and inflow from the aquifer of 5.5 mm. The average water balance of Lake Volvi was estimated at −38.89 mm/year.

The trophic status and the water quality of Lake Volvi are mainly affected by agricultural runoff, animal husbandry effluents, untreated or semi-treated domestic effluents and industrial wastewaters, mainly from food, dairy and other industries [52], as well as from suspended stream sediments and eroded bank materials.

Physical, chemical and biological variables for Lake Volvi were randomly measured for some short periods [52,53,54,55].

The mean values of SRP, TP and Chlα in water samples collected just below the surface of the water, as presented by Kaiserli et al. (2002) [52], were 60, 145, and 7.6 mg m⁻³, respectively. Similar values of SRP, TP and DO were presented by Petaloti et al. (2004) [55].

In a report of the Hellenic Ministry of Rural Development and Foods [56], a comparison of mean values of DO, phosphorus and nitrogen concentrations measured near lake surface water and near the bottom during two periods (1999–2000 and 2010–2012) was presented.

Table 2. Average values of some water quality parameters for the periods of 1999–2000 and 2010–2012.

	1999–2000	2010–2012	1999–2000	2010–2012
	Surface Water		Bottom Water
DO, mg L−1	8.2 to 9.8	8.5 to 10.4	6.0 to 6.4	2.1 to 5.1
P2O5, mg L−1	0.15 to 0.19	0.13 to 0.15	0.19 to 0.29	0.08 to 0.15
Inor N, mg L−1	0.10 to 0.98	0.11 to 0.39	0.17	0.49

shows the average values of these parameters.

More systematic measurements of temperature, dissolved oxygen and phytoplankton (as chlorophyll-α), and nutrient concentration were conducted under the National Water Monitoring Network program during the years 2011 to 2015 [37].

Thus, the calibration and recalibration of the model had to be based on the available data from 2013 to 2015. Simulations started on 1 January, when, based on measurements and meteorological data, the initial condition for water temperature was set at 9 °C. The initial condition for chlorophyll-α, SRP, OP and DO was defined from sporadic data from other years [56] and the fact that the lake is fully mixed during winter.

The values of the constants and coefficients in most simulations were based on a combination of the values given in the literature and on a sensitivity analysis. In this paper, the values taken from the literature were corrected during the calibration process [15,22,39,57,58,59,60,61,62,63].

3. Results and Discussion

The mathematical model QUALAKE was used to simulate the temperature and the concentrations of SRP, OP, Chlα, and DO in Lake Volvi on a daily basis for the years 2013 to 2015. Simulations were started on 1 January and continued in daily time steps until the end of each year.

In the following paragraphs the results of the temperature, DO, SRP, OP, and Chlα simulations will be presented and analyzed.

The calibration of QUALAKE was based on the available data from the years 2013 to 2015. During the calibration, by the method of trial and error, the model parameters were adjusted so that the simulated results of T_w, DO, SRP, OP and Chlα compared well with the measured values. The parameters that were used in the simulations are given in

Table 3. Values of constants and coefficients used in the simulation.

Symbol	Constants and Coefficients	Value	Units	Equation
βs	fraction of net short-wave radiation absorbed at the water surface	0.35		(8)
η	extinction coefficient for solar radiation	0.60	1/m	(8)
σ	constant in relation to turbulent diffusion coefficient and Richardson number	0.012		(8)
n	constant as σ	0.95		(9)
c2	constant in relation to reference value of turbulent diffusion coefficient and friction velocity	0.08		(9)
Phosphorus
apc	phosphorus to chlorophyll-α concentration ratio	0.5	mg P mg−1 Chlα	((11),(12))
kmin	mineralization constant	0.01–0.02	day−1	((11),(12))
kmchl	Michaelis–Menten constant for the mineralization of OP	25.0	mg Chlα m−3	((11),(12))
kSDp	release rate of SRP from the sediments	0.5	mg SP m−2 day−1
KDOSD	Michaelis–Menten constant for the release of SRP from the sediments in relation to the DO concentration	0.05	g DO m−3
wsop	vertical velocity of organic phosphorus sedimentation	0.005	m day−1	(12)
Phytoplankton
PHmax	maximum growth rate of phytoplankton	1.5	day−1	((11),(13))
qopt	optimal light intensity level	2.95	MJ m−2 day−1	((11),(13))
Kp	Michaelis–Menten constants for phosphorus	4	mg P m−3	((11),(13))
kr	coefficient of respiration loss of phytoplankton	0.07	day−1	((11),(13))
kmtot	coefficient of mortality loss (including grazing) of phytoplankton	0.05	day−1	(13)
ws	vertical velocity of phytoplankton sedimentation (sinking velocity)	0.005	m day−1	(13)

.

3.1. Water Temperature Simulation

In Figure 2, some characteristic simulated vertical lake water temperature profiles for specific days of the three years of the simulation are presented. Lake Volvi has the characteristic annual temperature cycle of temperate monomictic lakes. The model simulates well the onset of stratification, mixed layer depth and water temperature. In Figure 3, a comparison between simulated water surface temperature values and the air temperature is presented, while in Figure 4, a comparison of the simulated water temperature with the measurements at 2 m depth is presented.

The RMSE between the calculated and measured temperature values estimated for the available measurements of the three years of the simulation was 2.58 °C, while the coefficient of determination (r²) was 0.923. The RMSE value of temperature simulations for Lake Vegoritis (at a northern latitude of 40°47′) ranged from 0.42 to 1.84 °C [20]. The average RMSE and r² values of the water temperature simulation in nine lakes of different morphological characteristics were 1.1 °C (ranging from 0.68 to 1.93 °C) and 0.92 (ranging from 0.89 to 0.95), respectively [64].

According to simulation results, the lake showed a characteristic annual temperature cycle. Temperature stratification was observed from the beginning of April, with a thermocline depth of 3 to 8 m. The depth of the thermocline increased during the summer and remained deep until October. Fall overturn occurred between the end of October and the beginning of November. The maximum depth of the thermocline reached 9 m. The thermocline depth increased faster and was deeper for the year 2014, maybe due to the higher values of the diffusion coefficient as a result of the higher wind velocity of this year.

The model evaluated evaporation and all energy budget components of Lake Volvi on a daily basis, as well as the daily distribution of lake water temperature with depth, the surface temperature (T_sw), and the change in the thermal content of the lake (Q_t). The water temperature near the water–atmosphere interface of the lake controls the partitioning of available energy into latent and sensible heat [32]. An increased interest has developed in the water temperature and heat fluxes at the water–atmosphere interface during the last decades due to climate change. The results of energy budget and heat transport in the lakes, using models such as QUALAKE, could be considered reference data to develop different models based on artificial neural networks and energy budget models to estimate lake water surface temperature.

In Figure 5, the evaporation (E) values and thermal content change (Q_t) during the three successive years of simulation are presented. The energy budget components of Lake Volvi for each year of the simulation and their average values are given in

Table 4. Energy budget components (in MJ m⁻²d⁻¹) of Lake Volvi for the three years of simulation (2013 to 2015).

		Rs	Ra	Rbr	LE	H	Qt	E, mm day−1
2013	Ave	15.11	28.64	33.82	7.53	0.79	−0.32	3.07
	Max	28.39	36.11	39.39	23.90	6.05	11.11	9.73
	Min	2.62	12.27	28.70	−3.23	−10.17	−17.54	−1.31
2014	Ave	14.26	28.94	33.62	7.16	0.53	0.03	2.92
	Max	29.20	35.60	38.75	27.88	9.54	8.79	11.36
	Min	1.86	21.15	28.49	−7.51	−8.02	−23.34	−3.06
2015	Ave	15.26	27.88	33.10	7.60	0.76	−0.21	3.10
	Max	27.22	35.17	38.88	24.15	7.49	12.59	9.84
	Min	2.55	19.10	27.31	−5.13	−8.71	−17.99	−2.09
2013–2015	Ave	14.88	28.49	33.51	7.43	0.69	−0.17	3.03
	Max	29.20	36.11	39.39	27.88	9.54	12.59	11.36
	Min	1.86	12.27	27.31	−7.51	−10.17	−23.34	−3.06

Notation: R_s is the short-wave radiation; R_a is the long-wave radiation from the atmosphere; R_br is the back (long-wave) radiation; Q_t is the change in the energy content of the water body; H is the sensible heat; LE is the latent heat flux; and E is the evaporation (mm day⁻¹).

.

The average, maximum and minimum daily evaporation values are 3.03, 11.36 and −3.06 mm d⁻¹, respectively, and the thermal content change values are −0.17, 12.59 and −23.34 MJ m⁻² d⁻¹, respectively. The average evaporation was higher in 2013 and lower in 2014. The mean evaporation rate estimated, with the QUALAKE model, for Lake Vegoritis was 2.59 mm day⁻¹ [21].

Evaporation rates are low from the onset of calculations during winter and early spring. In summer, evaporation increases rapidly and remains high until early autumn. Negative evaporation rates (condensation) emerge during winter months.

3.2. Phytoplankton in the Lake

In Figure 6, the simulated values of the chlorophyll-α concentration during the three successive years of the simulation are presented.

The annual variation of chlorophyll-α concentration at four different depths, along with the measured values at the 2 m depth, is presented in Figure 7. A comparison of the simulated results from the QUALAKE model against the measured values of chlorophyll-α at a depth of 2 m shows generally acceptable results. The variation of chlorophyll-α concentration in the epilimnion (depths 2 and 5 m) presents different patterns along the different years. The concentration at 2 m depth shows higher variation during the years 2014 and 2015. There are some small periods in the spring months of 2015 with peaks in concentration. The distribution at a depth of 5 m was smoother throughout the study period. The RMSE between the calculated and measured values of chlorophyll-α estimated for the available measured data of the three years of the simulation was 10.26 mg m⁻³.

During thermal stratification, there are intense differences in the vertical distribution of Chlα concentrations regarding the epilimnetic and hypolimnetic values, as well as the depth of the maximum epilimnetic values. This depth depends on the depth of the epilimnion and the euphotic zone and the circulation and stratification conditions of the lake. The concentrations in the hypolimnion (depths 10 and 15 m) were, as expected, significantly lower in comparison to the values in the epilimnion and were zero after April.

The average Chlα concentration during the three years of the simulation was 4.44 mg m⁻³, with a minimum annual average of 1.31 mg m⁻³ in 2013 and a maximum of 11.28 mg m⁻³ in 2015 (

Table 5. Average, minimum and maximum values of water quality parameters in water surface and at 2 m depth and the average with depth values of SRP and Chlα from 2013 to 2015.

	Ta, °C	OP (0 m), mg m−3	SRP (0 m), mg m−3	Chlα (0 m), mg m−3	Tw (2m), °C	DO (2 m), mg L−1	Average SRP, mg m−3	Average Chlα, mg m−3
Ave2013	15.48	3.59	1.05	12.36	16.76	8.90	0.52	3.39
Max2013	29.30	4.27	2.27	31.46	27.51	12.29	0.71	7.86
Min2013	−2.50	2.32	0.44	5.85	6.47	6.67	0.42	1.31
Ave2014	15.68	5.01	1.13	17.96	16.44	9.12	0.54	4.22
Max2014	29.10	5.75	3.18	40.82	27.04	12.32	0.79	8.70
Min2014	−2.50	3.62	0.44	6.15	6.74	6.75	0.44	1.63
Ave2015	14.24	5.72	1.53	21.92	14.75	10.57	0.67	5.71
Max2015	28.80	6.27	6.58	42.41	25.34	13.18	0.86	11.28
Min2015	−4.10	5.19	0.47	8.67	5.20	7.34	0.57	3.16
ave2013–2015	15.13	4.77	1.23	17.41	15.99	9.53	0.58	4.44
max2013–2015	29.30	6.27	6.58	42.41	27.51	13.18	0.86	11.28
min2013–2015	−4.10	2.32	0.44	5.85	5.20	6.67	0.42	1.31

). In general, the Chlα concentrations were higher in 2015.

3.3. Phosphorus in the Lake

There are no measured data of phosphorus in Lake Volvi. Some measurements of total phosphorus concentration have been conducted on behalf of the Hellenic Ministry of Rural Development and Foods [56] and the National Water Monitoring Networks [37] for some specific years, but only in the surface layer of the lake.

A comparison of the SRP profiles simulated on different days of the three years is presented in Figure 8. In the spring, the distribution was similar under the thermocline depth, though very different in the surface layer. In the 3 to 8 m layer (metalimnetic layer), a decrease in concentration was observed during July, mainly due to the high phytoplankton production in this layer. Later in October, significantly lower surface concentrations were estimated in comparison to the previous months, as phytoplankton production remains high during this period of the year. The thermocline inhibits the mixing of the two layers (epilimnion and hypolimnion), and, therefore, SRP, which is produced in the hypolimnion by the mineralization of the sinking dead organic matter that accumulates there and continuously increases hypolimnetic SRP concentrations. In the deeper layers, the SRP concentrations for the years 2014 and 2015 were similar, although they were lower in 2013, mainly due to the lower concentrations of Chlα for this year.

The temporal distribution of SRP concentration during the three years is shown in Figure 9. Very intensive variation was observed at the surface, with high peaks in 2014 and especially in the first months of 2015. There were some periods in which the concentration of SRP was very high. These are connected to high rainfall and high runoff, which enriched the lake with inorganic phosphorus. Deeper in the hypolimnion, the calculated concentrations of SRP were relatively constant and, in general, lower than 2 mg m⁻³, although at different sizes during the different years of the simulation. In the 3 to 5 m depth layer near the thermocline, the SRP concentrations presented lower values in relation to the concentrations at other depths, as this is the layer where maximum phytoplankton growth takes place.

The cumulative daily loading of SRP in the lake for each year and the daily rainfall of the three-year period are presented in Figure 10. The total SPR loading of the three years was 85.65, 128.31 and 111.65 mg m⁻², respectively, while the cumulative rainfall was 420.6, 914.7 and 573.6 mm, respectively.

The average SRP concentration in the lake during the three years of the simulation is presented in Figure 11. The mean value of the three years was 0.57 mg m⁻³, with a minimum average of 0.40 mg m⁻³ in 2013 and a maximum of 0.85 mg m⁻³ in 2014 (Table 5). A comparison of Chlα and SRP average concentrations (Figure 11) shows that there is an inverse correlation between the values of the two parameters.

The annual variation of OP concentrations at specific depths is presented in Figure 12. The concentration of OP in the epilimnion presented significant variation throughout the year. There is an increase in OP concentrations after February until the end of May, followed by a decrease with the lower values at the end of July. Then from September onward, an increase was again observed. The OP concentrations ranged from 2.3 to 4.2 mg m⁻³ in 2013, were higher in 2014 from 3.6 to 5.7 mg m⁻³, and from 5 to 6.3 mg m⁻³ in 2015 in the epilimnion and lower but continuously increased year to year in the hypolimnion (2.5 to 3.8 in 2013, 3.1 to 5.5 in 2014, and 5 to 6 mg m⁻³ in 2015).

Throughout the year, near-bottom OP concentrations were higher in relation to hypolimnetic concentrations. During the period of thermal stratification, OP is rapidly recycled in the epilimnion because of the high concentrations of phytoplankton and, as a result, zooplankton organisms. High OP values are observed in the bottom of the hypolimnion and the sediments, where the particulate portion of OP sinks and decomposes relatively slowly. The OP concentrations near the bottom throughout the year are higher in relation to hypolimnetic concentrations. The hypolimnetic OP concentrations were similar in 2013 and 2014, while they were notably higher in 2015, mainly due to the high phytoplankton growth rate and Chlα concentrations of this year.

The cumulative loading of OP to the lake is presented in Figure 10. The total OP loading for the three years of 2013–2015 was 12.36, 18.52 and 16.11 mg m⁻², respectively.

3.4. Dissolved Oxygen Simulation

In Figure 13, some characteristic simulated DO profiles are presented for specific days from the three years. DO concentration profiles during the first months of the year were nearly homogeneous, while during summer and early fall showed different distributions along the depth for each year of the simulation, following the patterns of the Chlα and temperature distributions. In the euphotic zone, the DO concentration decreased throughout each year. The DO values were less than 2 mg L⁻¹ in the hypolimnion during the autumn of 2014, while the DO concentrations in 2013 and 2015 remained higher than 4 mg L⁻¹. Due to stratification, the hypolimnion is cut off from the high epilimnetic oxygen concentrations. Non-occurrence of the metalimnetic oxygen minimum and the DO concentrations in the hypolimnion are evidence that lake Volvi is not under eutrophication pressure.

In Figure 14, the daily DO concentration at 2 m depth in comparison to the measurements taken during the three successive years is presented. The RMSE between the calculated and measured DO values estimated for the available measurements of the three years of the simulation was 2.085 mg L⁻¹. The average RMSE values of DO in a number of lakes in the USA, presented by Stefan et al. (1993), were from 0.6 to 2.3 mg L⁻¹ with an average of 1.4 mg L⁻¹ [64].

3.5. Discussion

The evaluation of the model was based on measured data of temperature, DO and Chlα during three years from 2013 to 2015 in lake Volvi. This time period was selected as a result of the available data from the National Observatory of Athens in Greece [37]. Assessing the model based on the values of r² and RMSE for these parameters, the performance of the QUALAKE model in Lake Volvi was quite good, especially when taking into consideration the uncertainties in model parameters and nutrient loading estimations and the absence of measurements of other water quality parameters. The values of the performance criteria are generally similar to those reported in other ecosystem model studies [10,64,65,66].

Parameters estimation is a critical point in the development of environmental models, according to Chapra (1997) [42]. The identification of parameter values is a very difficult task due to the inherent complexity of environmental systems, scarcity of data, and mathematical and computational limitations [67]. Saloranta and Andersen (2007) [8] noted that the performance of a model code is very much dependent on the particular task to which it is applied, implying that a general model code “validation” may be an unrealistic goal based on graphical estimates.

The energy budget for a lake relates to the net transfer of energy to and from the water. The temperature of the water and its variations express the changes in the stored thermal energy, and it is one of the principal water quality and quantity factors of lake hydro ecosystems. The changes in lake water temperature and the temperature stratification dynamics can have a profound effect on lake biological and chemical processes, as well as on their hydrological cycle.

Incoming (R_a) and outgoing (R_br) long-wave radiation are the energy budget components with the largest mean annual rates (R_a = 28.49 MJ m⁻² day⁻¹ and R_br = 33.51 MJ m⁻² day⁻¹). Because of their opposite signs, there is a mean loss of energy of 5.12 MJ m⁻² day⁻¹ (Table 4). Therefore, the source of energy is the net short-wave solar radiation (R_s) with a mean rate of 14.88 MJ m⁻² day⁻¹. R_s, R_a and R_br radiation values exhibit a characteristic seasonal variation, with the highest rates in summer and the lowest in winter. Similarly, net radiation (R_n) obtains its largest values in the summer and its lowest in the winter, with a three-year average value of 7.97 MJ m⁻² day⁻¹. Its mean annual value is explicitly positive, as net radiation is the source of energy for all the physical and biological activities of the lake ecosystem.

The seasonal variation of the thermal content of the water body (Q_t) is balanced between its positive and negative values of incoming and outgoing radiation and sensible and latent heat (Figure 15). In spring, the lake consecutively gains and stores thermal energy through net radiation. In summer, the energy budget is positive. In autumn and winter, the energy budget is negative since the stored energy is released as sensible and mainly as latent heat, resulting in the gradual cooling of the lake. For the three years of simulation, the total energy gain of Lake Volvi was −0.17 MJ m⁻² day⁻¹.

The average SRP concentrations were low for all years. SRP is rapidly used for phytoplankton growth in the euphotic zone due to the high primary productivity of the lake. The sinking dead organic matter accumulates in the hypolimnion, and its mineralization increases SRP concentrations in this layer, which cannot be transferred to the epilimnion because of the thermocline. The different meteorological conditions, with higher rainfalls at the end of 2014 and the beginning of 2015, had an important effect on inorganic phosphorus loads, increasing the SRP concentration in the lake and consequently increasing the photosynthesis rate. The latter also affected the DO concentrations in the epilimnion during the same period.

The average annual OP concentrations in the hypolimnion ranged from 2.8–5.4 mg m⁻³ during the three years. Similar, but slightly higher, were the values of OP in the epilimnion, which ranged from 3.6 to 5.7 mg m⁻³. The distribution with depth of OP values was nearly the same for all the years of the simulation, with only higher near the bottom values during 2015. A specific portion of OP sinks into the hypolimnion and the sediments, where it decomposes relatively slowly.

The Chlα concentration differed from year to year in regard to both magnitude and pattern (Figure 7). The maximum average Chlα concentration in the lake ranged from 7.79 mg m⁻³ (2013) to 10 mg m⁻³ (2015). The peaks in the Chlα concentration appeared from mid-spring to mid-summer, while in 2015, three peaks were observed (spring, summer, and autumn), which is characteristic of lakes that are becoming more productive [68]. The average Chlα concentration agrees with the WFD’s (the EU Water Framework Directive) water quality requirements for good ecological quality (Chlα concentration of deep-water lakes should not exceed 12 mg m⁻³).

The main environmental variables affecting primary productivity are nutrient availability, temperature, and solar radiation. As phosphorus is considered the limiting nutrient for phytoplankton growth in most lakes, phytoplankton concentrations are usually closely related to SRP concentrations. Comparing the distribution of SRP and Chlα from Figure 7, Figure 9 and Figure 11, the results show that phytoplankton concentrations throughout the year are strongly affected by SRP concentrations and vice versa [22]. Thus, high primary productivity at the end of winter and early spring causes concentrations of SRP to decrease rapidly during spring, followed by a reduction in phytoplankton production rates. During summer, a temporary recovery in SRP is observed as a result of a higher rate of mineralization due to the presence of dense populations of zooplankton organisms, which accelerate the process of mineralization. Huang et al. (2012) [69] considered that the chlorophyll-α concentration was highly sensitive to parameters directly related to growth rates (such as maximum phytoplankton growth rate and temperature, nutrient, and light factors).

The model results of the DO simulation showed that there was good agreement between the simulated and measured values of DO for the available data at a 2 m depth in the lake. According to measurements from the Hellenic Ministry of Rural Development and Foods [56], the mean values of DO concentrations measured near the lake water surface ranged from 8.2 to 9.8 during 1999–2000 and 8.5 to 10.4 during 2010–2012, and near the bottom of the lake from 6.0 to 6.4 during 1999–2000 and 2.1 to 5.1 during 2010–2012 (Table 2). Figure 14 and Figure 15 showed a similar distribution of DO in Lake Volvi.

Suitable calibrated and verified models for energy and water quality of lake ecosystems are necessary for lake water and ecosystem management. A lake’s ecosystem responds directly to climate change. Increasing water temperature affects the hydrodynamics in lakes, expanding the thermal stratification period and the depth of the thermocline [65,70,71]. Water temperature, the most sensitive parameter to climate change, is a key parameter in most biological reactions, influencing water chemistry, biochemical reactions, and aquatic organisms directly.

The results of energy and water quality parameters estimated for Lake Volvi for the three years of simulation showed that the model could be used for the management of the water quality of the lake. It can also be used to forecast the changes and effects in the status of the lake under climate change scenarios. The model should be recalibrated with more recent measured data either from the National Water Monitoring Network or more intensive measurements in the lake.

4. Conclusions

The one-dimensional mathematical model QUALAKE was used to describe energy (temperature, heat balance and evaporation), water quality (dissolved oxygen, inorganic and organic phosphorus) and biological productivity (phytoplankton–chlorophyll-α) during the annual cycle in a stratified lake and the results were presented and evaluated. The model simultaneously describes water temperature and concentrations of DO, SRP, OP, and phytoplankton (as chlorophyll-α) in a daily time step. It simultaneously calculates evaporation rates and energy budget components.

Data from Lake Volvi in Central Macedonia of Greece (in northern Greece) for three successive years (2013 to 2015) were used to calibrate and recalibrate the model.

The simulation results of the model’s application are an important source of information on changes in the physical, chemical and biological variables of the lake for which very few data are available. The energy, nutrients, and biological submodels were parameterized using literature data and verified by computing simulations of temperature, DO and Chlα patterns with field data.

The model results of water temperature, dissolved oxygen and chlorophyll-α were compared to measurements at 2 m depth for the years 2013 to 2015. The comparison showed that the predicted values of these parameters were all in good agreement with the measurements.

The results of the model simulations during the three successive years for energy, water quality and primary productivity showed that the average concentration of the studied parameters was 0.58 mg m⁻³ for SRP, 4.5 mg m⁻³ for OP, and 4.4 mg m⁻³ for Chlα. The DO in the epilimnion was near saturation, while in the hypolimnion remained higher than 2 mg L⁻¹. The average evaporation from the water surface was 3.05 mm day⁻¹.

There is a general deficiency in water quality modeling in Greece, mainly due to the lack of appropriate data for the application and validation of such models. This is the first time that a mathematical model, such as the QUALAKE model, has been applied to Lake Volvi, one of the most important impoundments in northern Greece, protected by the Ramsar International Convention. Furthermore, according to the EU Water Framework Directive 2000/60, Annex II 1.5, EU member states: “shall carry out an assessment of the likelihood that surface water bodies within the river basin district will fail to meet the environmental quality objectives set for the bodies under Article 4. Member states may utilize modeling techniques to assist in such an assessment”. The QUALAKE model could be used in fulfilling the Directive requirements, regarding Lake Volvi, in view of the scarcity of water quality data for this lake and the fact that a management plan for the river basin district of Lake Volvi has not yet been elaborated.