Mathematical modelling as a tool for the assessment of impact of thermodynamics on the algal growth in dam reservoirs – case study of the Goczalkowice Reservoir

Changes in water temperature in lakes and reservoirs is one of the key processes affecting the spatial and temporal distribution of chemical compounds and biological elements such as plankton and fish [Elci 2008; Genova et al. 2010; Yu et al. 2010]. The water temperature depends on geographic conditions (e.g., bathymetry, elevation, latitude and interactions with groundwater) and seasonal changes in meteorological conditions (e.g., solar radiation, cloud cover, wind speed and direction) [Lee et al. 2013; Tuan et al. 2009; Tarasiuk et al. 2015]. The temperaturebiochemistry interdependences are often bidirectional, and the concentration of suspended solids or plankton inhibiting the penetration of solar radiation, for example, may impact the water temperature [Kumagai et al. 2000; Mazumder et al. 1990]. Temperature affects the chemical and biological status directly, as particular processes occur more intensively at a given range of temperature and indirectly influence the water density, and thus, water mixing or stratification [Elci 2008]. In addition to the effects of temperature, water mixing processes also depend on the inflow and outflow currents, wind and rainfall. The strength of all these factors depends regionally on the climate and geographical conditions (e.g., rainfalls and currents can dominate in monsoon seasons and inflow currents can be the main factor during heavy rainfalls in mountains or foothills areas) [Heungsoo et al. 2013].

During summers, the warm water of epilimnion (upper well-mixed zone) is unable to mix with more dense water in hypolimnion (lower, colder zone) because of density differences that cannot be disturbed by the energy of winds, inflow or outflow currents and rainfall [Boehrer, Schultze 2008; Elci 2008; Hodges 2000; Lee et al.2013; Yu et al. 2010]. When the density differences arise in water column, the upper layer is usually influenced by winds and therefore well mixed. In contrast, the lower zone is cut-off from atmospheric oxygen and light. Photosynthesis may be inhibited, and the oxygen consumed by bacteria and other organisms. Near the bottom, anoxic conditions may occur. In such conditions, a solubility of nutrients, heavy metals and hydrogen sulphides increases and poses a risk of release from bottom sediments [Boehrer, Schultze 2008; Hudson, Kirschner 1997; Lee et al. 2013]. In addition to release from sediments, the organic matter falling from the epilimnion increases nutrients accumulation [Yu et al. 2010]. During the fall turnover, water mixing status changes as the solar radiation energy is outweighed by winds, currents and other factors causing vertical water flows. During the more intensive water mixing, substances accumulated in deep (hypolimnietic) water become available near the surface [Boehrer, Schultze 2008]. In addition to the seasonal changes in reservoir thermodynamics, shallow ones can stratify and destratify on a daily basis [Branco, Torgersen 2009].

Analyses of hydrodynamics, thermodynamics and impacts on water quality and ecosystems can be supported by mathematical modelling. There is a variety of numerical models describing the physical processes affecting water mixing and heat transport or coupled hydrodynamic and water quality models. Jensen et al. (2015) reported over 1500 such tools based on a literature review. There are numerous one-dimensional models (e.g., CE-QUAL-R1, DYRESM, DUFLOW, GLM, GOTM, LIMNMOD, MINLAKE, Mylake, PROTECH, SIMSTRAT) widely applied and in most cases, allowing for analyses of vertical water mixing and/or quality [Jensen et al. 2015]. These models were often a basis for two dimensional models such as CE-QUAL-W2, which can be used to simulate changes along rivers and narrow reservoirs, however, with the assumption that horizontal variations in water temperature and quality are of minor importance (not valid for analyses of local algal blooms) [Romero et al. 2004; Lee et al. 2013]. The most advanced tools are three dimensional models (e.g., AEM3D, CAEDYM, ELCOM, GEMSS, GETM) or modelling systems that allow to choose the number of dimensions according to the needs of application (e.g., Delft3D, EFDC, WASP) [Jansen et al. 2015]. Analyses of advantages and disadvantages of various models and modelling approaches were a subject of numerous studies (e.g., Saloranta et al. 2004; Gao, Li 2014; Jensen et al. 2015).

In the presented study, the ELCOM model was applied and calibrated (in terms of the water temperature) in order to simulate the summer and fall changes in water temperature of the Goczałkowice Reservoir and to prepare a basis for the comparisons of the reservoir’s thermodynamics and its chemical and ecological status. Mathematical model used and described in this study was applied as a part of task 5: ‘Development, application and verification of models for lakes physics, chemistry and ecosystems’ of the ZiZOZAP project.

Goczałkowice Reservoir is the biggest dam reservoir in the south of Poland. It covers over 32 km² and has two main inflows: the Small Vistula River (the main inflow, approx. 80%) and the Bajerka River with total basin area of 530 km² (Figure 1).

Figure 1

The reservoir was created in 1955 and serves as a major part of the system supplying the Upper Silesian agglomeration (approx. 3.4 million inhabitants) with potable water. It is also a storage reservoir protecting downstream areas from floods and droughts. Additionally, the reservoir being a part of the Natura 2000 system helps to protect a wide range of habitats and species [Dabioch et al. 2013; Młynarczyk et al. 2013; Polak et al. 2011], The reservoir is included in the national monitoring system and since 2010, extensive research monitoring has been carried out in the framework of the ZiZOZap project [ZiZOZap 2010], Both the operational monitoring system and real-time research measurements indicate strong changes in the water quality posing a risk to the water treatment plants and biodiversity dynamics (parameters of concern include: nutrients, sediments, dissolved oxygen and heavy metals). Some of the water quality parameters (e.g., phosphorus, humic acids and heavy metals) were recently analysed in relation to their abilities to be released from sediments [Dabioch et al. 2013; Młynarczyk et al. 2013; Polak et al. 2011]. However, the reservoir hydrodynamics and thermodynamics were not analysed in order to identify their impact on concentrations of chlorophyll a. To assess the interactions between reservoir thermodynamics and the concentration of chlorophyll a, the first step was to obtain detailed information on the changes in water velocity and temperature. For this purpose, a mathematical model of the reservoir was applied.

In this study, the Estuary, Lake and Coastal Ocean Model (ELCOM) was applied. ELCOM is a three-dimensional hydrodynamics model for lakes and reservoirs and is used to predict the variation of water temperature and salinity in space and time. Heat exchange through the water’s surface is governed by standard bulk transfer models. The energy transfer across the free surface is separated into non-penetrative components of long-wave radiation, sensible heat transfer, and evaporative heat loss, complemented by penetrative shortwave radiation. Non-penetrative effects are introduced as sources of temperature in the surfacemixed layer, whereas penetrative effects are introduced as source terms in one or more grid layers on the basis of an exponential decay and an extinction coefficient [Hodges, Dallimore 2013].

The model of Goczałkowice Reservoir presented here was applied (tested and calibrated) for a period of 6 months: June – November 2010. The time step used in simulations was set to 5 minutes. The starting time was determined by the setup of the real-time monitoring system, which was fully operational in June. The six-month period of analyses was chosen in order to be: 1) long enough to include processes crucial from the ecological and management points of view (e.g., algal bloom, high concentrations of suspended solids) and 2) as short as possible to enable multiple iterations required for the model calibration in an acceptable time.

The model of Goczałkowice Reservoir consists of 20 layers of thickness varying from 0.5 to 1.25 m. The thickness of 0.5 was applied to the surface water layers (up from 255.25 m a.s.l.), the thickness of 1.25 m was applied to the bottom layer (244–245.25 m a.s.l.) and the thickness of 1 m was used for ten intermediate layers. Twelve of these 20 model layers are usually wet, with the average water table elevation at 255.5 m a.s.l. The next three layers of 0.5 m thickness are reserved for higher water levels, with the maximum impoundment level at 257 m a.s.l. The top five layers should be always dry and were left only for calibration purposes. Horizontal resolutions of 25, 50, 100 and 250 m were tested and finally a resolution of 100 m was chosen as a compromise on the level of details and the calculation time (one ELCOM simulation of 180 days lasted 35 minutes for 250 m resolution, 5 hours for 100 m resolution and 9 days for the resolution of 25 m). The chosen resolution resulted in 119 columns, 60 rows and 41,683 total calculation cells. The model includes 7 inflows: the Vistula River (the main inflow), the Bajerka River and five pumping stations transferring the water excess from depressed agricultural and grassland areas around the western part of the reservoir (Figure 1). The average inflows in the analysed period are as follows: Vistula 7.98 m³/s, Bajerka 0.40 m³/s and pumping stations 0.41 m³/s. Outflows include: intake (average 2.02 m³s), spillway (average 7.07 m³/s) and bottom outflow required for the protection of ecosystems downstream (0.6 m³s). The temporal resolution of all inflow and outflow data is one day and all inflows include the information on water temperature.

The model included inputs from one meteorological station located above the water surface near the Bajerka river inflow. Meteorological data are of hourly resolution and are presented in Figure 2. The real-time water monitoring system launched in 2010 includes hourly measurements of water temperature in one water column with one-meter intervals (Figure 2). These measurements were the main input for the model calibration and verification. Additionally, the water temperature was measured periodically (approx. once a month) in 7 research monitoring points and in two monitoring points of the national monitoring system (Figure 3). These data were used only for model verification purposes.

Figure 2

Figure 3

The surface heat exchange includes shortwave radiation, longwave radiation, sensible heat flux and evaporation. The latter three are constituents of the non-penetrative energy density deposited in the surface layer. Solar radiation is one of the most important factors affecting the lake hydrodynamics and ecosystem. The ELCOM model divides the shortwave radiation into four bands: Photosynthetically Active Radiation (PAR), Near Infrared (NIR) and two Ultra Violet bands: A (UVA) and B (UVB). Percentages for each band in the whole shortwave radiation were set in the model as follows: PAR = 45%, NIR = 51%, UVA = 3.5% and UVB = 0.5%. The depth of the shortwave penetration depends in the model on the extinction coefficient, which is set separately for each band (Table 1).

Table 1

Parameters used in ELCOM model.

The penetration was calculated according to the Beer-Lambert law [Hodges, Dallimore 2013]. In case of an excess of shortwave energy at the bottom of water columns, 90% of the energy is allowed to propagate back. For the longwave energy, two options were tested for the Goczalkowice Reservoir model: first based on the albedo and longwave radiation energy density deposited into the surface layer of water, and second, based on the albedo, cloud cover fraction and the temperature of the air and surface water layer [Hodges, Dallimore 2013]. The sensible heat loss from the surface of the reservoir was calculated based on the wind speed (at 10 m reference height above the water surface), the difference between the surface water temperature and air temperature, and the sensible heat transfer coefficient (Table 1). Finally, the last component of non-penetrative energy, the evaporation, was calculated based on atmospheric pressure and density, vapour pressure, saturated vapour pressure, water surface temperature, wind speed and heat transfer coefficient. The variability in the sensible and latent (evaporation) heat transfer described above can be affected by air column stability and water roughness [Imberger, Patterson 1990]. These factors have the effect of altering heat exchange coefficients, especially if meteorological parameters are measured over the water surface, as it is in the case of Goczalkowice Reservoir. ELCOM includes an atmospheric stability correction procedure applied after Hicks [1975] [Hodges, Dallimore 2013].

The surface thermodynamics described above are applied in ELCOM to the upper grid cells changing their temperature. Heat flux is added to the uppermost grid cell while solar radiation is added to the water column using exponential decay over depth. The heat transfer changes the density stratification, which is also affected by the concentration of suspended solids. Once the new density field is calculated, the mixing process is modelled on a layer-by-layer basis through each water column [Hodges, Dallimore 2013].

The model of Goczalkowice Reservoir excludes: the Coriolis term, reservoir generated tides and rainfall input. Whilst the excluding of the Coriolis term and tides may be easily explained by their minor importance for the analysed problem, the same could not be said of the rainfall input. Reasons for ignoring the rainfall effect on the reservoir thermodynamics are: 1) the geographical characteristics of the reservoir catchment (in this foothill area, the inflow currents and winds are the most dominant factors affecting the water mixing process), and 2) the observed low intensity of rainfalls in the analysed period. More about the factors determining the reservoir hydrodynamics and thermodynamics can be found in the ‘Results and discussion’ section.

ELCOM (and its successor – AEM3D) is a widely used model with teens of scientific papers confirming its applicability worldwide. Numerous applications are reported in Europe [Papadimitrakis, Karalis 2009; Lang et al. 2010; Imberger et al. 2017], and several in Poland, where it was used mainly to analyse the shallow dam reservoirs [Kliś et al. 2014; Karpiński, Łozowski 2017; Woźnica et al. 2017; 2017b]. Calibration of the Goczałkowice Reservoir model included: 1) modification of heat transfer coefficients, light extinction coefficients and initial water temperature in the reservoir or inflows, 2) using different meteorological data sources and 3) using different methods of longwave energy input to the surface layer (as described in the section above). The calibration process included over 400 iterations and the simulated water temperatures were compared to the output from the real-time monitoring system (Figure 4). It is a typical approach to calibrate the model’s thermodynamics prior to the application of the model’s results to the analyses of water quality or ecosystems [Zhang et al. 2008]. Such calibration and subsequent validation are often (but not as a rule, e.g., Leon et al. 2006) concluded using a statistical assessment of the model performance [Abbasi et al. 2016]. In case of the Goczalkowice Reservoir, model validation was based on both real-time and periodical monitoring and was evaluated using the determination coefficient (R²) and the Nash-Sutcliffe (NS) model efficiency coefficient [Nash, Sutcliffe 1970]. Water profile temperatures for the simulation before and after the calibration are presented in Figure 5. Before the calibration, R² = 0.68 varying at different water depths from 0.4 to 0.86, while NS = 0.46 varying from -0.22 to 0.86. Both parameters indicate that the model before the calibration explained less than half of the temperature variations in the reservoir. Moreover, the negative NS suggests that using an average temperature at several water depths could give a more accurate result than the ELCOM model. The calibration process resulted in significant improvement of the model efficiency. For the best simulation, R² is equal to 0.96, varying at different water depths from 0.9 to 0.96, while NS=0.95 varying from 0.92 to 0.97.

Figure 4

Figure 5

From the beginning of summer (middle of June) to late summer (end of August), according to the monitoring data, the average water temperature in the surface layer is approximately 3.3°C higher than at the bottom of the reservoir, whereas the maximum difference (ΔT) reaches 10°C (Figure 2). Water temperatures simulated by ELCOM are almost the same in the surface and bottom layers near natural banks and in the shallow parts of the reservoir (south and west). Temperatures differ in the deeper, central part of the reservoir and in the vicinity of former river channels. Despite this, the reservoir is considered as shallow and stratified in winters only [Czaplicka-Kotas et al. 2012], the difference of temperature near surface and bottom causes (in the model) differences in the water density, thus impeding the vertical flow. During the summer, three water mixing events occurred in the reservoir that were confirmed by almost constant temperature in water profiles (0–2°C differences) and by increased vertical water velocity calculated by model (example on Figure 6). The first event started on 24^th June and on 28^th June the temperature of water was uniform. During those days, the maximal difference in bottom and surface temperature (ΔT) decreased from 10 to less than 1°C in almost the whole reservoir. Similar events occurred in periods between 3^rd and 6^th August and between 14^th and 19^th August. These two events resulted in a decrease in ΔT from 4 to 1.5°C and from 5 to 0°C respectively. After the last mentioned short summer mixing events, a period considered as autumn occurred. At the end of August, variations in water temperature were less than 2°C. At the end of September, ΔT decreased further and remained near to 0°C for the rest of the analysed period.

Figure 6

The simulation of water temperature was followed by an analysis of relations between water temperature and observed concentration of chlorophyll a. As mentioned in the Introduction section (or as mentioned above), during summer and autumn, a strong variability in the concentration of chlorophyll a is observed in the Goczalkowice Reservoir. Apart from the period analysed in this study, it is also confirmed by other researchers, for example, Gwiazda et al. (2014), that concentrations range from 3.2 to 107.8 μg chla / l in 2011 and 2012. ChlA is widely considered as an indicator of algal biomass and of algal blooms, when high concentrations of the chlorophyll occur. Such blooms occur in the analysed reservoir with cyanobacteria, chlorophytes and diatoms as dominating groups [Czaplicka-Kotas 2012]. High concentrations of ChlA pose a risk to the water treatment process and decrease the cost-effectiveness of water supplies. Usually in lakes and reservoirs, a decrease in water temperature variations is observed when water is mixed. Such water temperature status allows the release of nutrients stored in the deeper layers of water as the water density becomes uniform. The result of increased concentrations of nutrients close to the water surface is an increase in phytoplankton growth rate. Comparisons made for the Goczalkowice Reservoir confirm that decreases in the temperature variations in the water profile are followed by increases in ChlA concentration. The analysis indicated that the average lag in ChlA response to the temperature in reservoir is about 60 hours. The moderate correlation of 0.44 (Pearson correlation coefficient – R) was observed across the whole analysed period, increasing to 0.6 in autumn. In summer, each of the three short events of water mixing was followed by an increase in the concentration of ChlA. The concentration of chlorophyll increased from 7 to 70 μg/l, from 7 to 90 μg/l and from 7 to 55 μg/l for the three events mentioned above (Figure 7). In autumn, water in the Goczalkowice reservoir is mixed, so further releases of nutrients from deeper layers of water do not occur. However, minor increases in ChlA concentration still occur along with increases in surface water temperature (result of short periods of increased air temperature and solar radiation).

Figure 7

The presented study confirms the applicability of the ELCOM model to simulate thermodynamics of a dammed reservoir using the Goczalkowice reservoir as an example. Simulated and measured water temperatures are matched with high accuracy (R² = 0.96, Nash-Sutcliffe coefficient = 0.95). The model allows assessment of changes in water temperature with high spatial and temporal resolution. In case of the Goczalkowice reservoir, it enabled identification of parts of the reservoir, where the temperature is more uniform or diversified in the water profile. The model also allowed to identify where and when vertical water mixing occurs. Such variations in temperature are difficult to be identified based on the monitoring system and even real-time measurements are spatially limited and should be complemented by modelling systems such as presented in this study.

Basing on the high correlation between water temperature (water mixing processes) and concentration of chlorophyll a, it is possible to identify areas prone to algal blooms with the same spatial and temporal resolution as used in the mathematical model of water thermodynamics. Moreover, the application of model confirmed that changes in ChlA concentrations are affected by the temperature of particular water layers and that each water mixing event leading to the unification of temperature in water profile was followed by a significant increase in algal growth (peaks of approx. 10 times the initial ChlA concentration). Such possibilities can be useful not only for the lake or reservoir characterisations but also for warning systems informing about the potential algal growth as a response to forecasted weather conditions.