3D-modelling of the thermal circumstances of a lake under artificial aeration

A 3D-model was developed to study the effects of hypolimnetic aeration on the temperature profile of a thermally stratified Lake Vesijärvi (southern Finland). Aeration was conducted by pumping epilimnetic water through the thermocline to the hypolimnion without breaking the thermal stratification. The model used time transient equation based on Navier–Stokes equation. The model was fitted to the vertical temperature distribution and environmental parameters (wind, air temperature, and solar radiation) before the onset of aeration, and the model was used to predict the vertical temperature distribution 3 and 15 days after the onset of aeration (1 August and 22 August). The difference between the modelled and observed temperature was on average 0.6 °C. The average percentage model error was 4.0% on 1 August and 3.7% on 22 August. In the epilimnion, model accuracy depended on the difference between the observed temperature and boundary conditions. In the hypolimnion, the model residual decreased with increasing depth. On 1 August, the model predicted a homogenous temperature profile in the hypolimnion, while the observed temperature decreased moderately from the thermocline to the bottom. This was because the effect of sediment was not included in the model. On 22 August, the modelled and observed temperatures near the bottom were identical demonstrating that the heat transfer by the aerator masked the effect of sediment and that exclusion of sediment heat from the model does not cause considerable error unless very short-term effects of aeration are studied. In all, the model successfully described the effects of the aerator on the lake’s temperature profile. The results confirmed the validity of the applied computational fluid dynamic in artificial aeration; based on the simulated results, the effect of aeration can be predicted.


Introduction
Hypolimnetic anoxia is a common phenomenon in productive lakes and reservoirs that are deep enough for thermal stratification. Anoxia can develop during the summer, because depletion of dissolved oxygen in the hypolimnion is much faster than transfer of dissolved oxygen (DO) from the surface layers by advection and dispersion (Mackenthun and Stefan 1998;Beutel et al. 2007). Anoxia accelerates eutrophication, because in low oxygen conditions, iron (Fe) is reduced and loses its capability to bind phosphorus (P) (e.g. Boström et al. 1982). This is of crucial importance for lake ecosystems, where P is one of the most important factors regulating ecosystem productivity and sediment has a large storage of P compared to the water column (Wetzel 2001). Due to the effect of anoxia on internal P loading, artificial aeration or oxygenation (for simplicity, the term aeration is used hereafter) is frequently used in the management of eutrophicated lakes and reservoirs suffering from oxygen deficits (Cooke et al. 2005;Kuha et al. 2016;Grochowska et al. 2017;Horppila et al. 2017). However, evidence on the long-term positive effects of aeration on water and sediment quality is obscure (Gächter and Wehrli 1998;Liboriussen et al. 2009;Kuha et al. 2016;Horppila et al. 2017). In addition, an important issue is that aeration has negative side-effects that can counteract the positive effects. These include increased temperature and water flow velocities in the hypolimnion with consequences for the oxygen demand of the sediment (Beutel et al. 2007;Gantzer et al. 2009;Gerling et al. 2014;Niemistö et al. 2016). Such effects complicate the planning of aeration projects and scaling the use of aerators. Numerical models for lake hydrodynamics have previously been used in studying the design of management methods. Three-dimensional models with complex turbulence modules have been developed for instance for the influences of hypolimnetic jets (Toffolon and Serafini 2013), simulating lake mixing by air-bubble plumes and evaporation effects (Helfer et al. 2011) and sediment transport (Hu et al. 2009). Still, no three-dimensional model has been presented for hypolimnetic aeration with units that pump epilimnetic oxygen-rich water through the thermocline to the hypolimnion without breaking the thermal stratification. Such systems have been, however, used in hundreds of lakes (Cooke et al. 2005). Computational fluid dynamic (CFD) is an accurate way of simulating water flow and water energy transport in engineering (Choi et al. 2013). In this paper, temperature was chosen as an indicator, and the CFD method was applied in modelling of artificial lake aeration. A 3-D model was set up and temperature distribution of an aerated stratified lake was calculated by this model. The accuracy of the model was certified by field data.

Study site
Enonselkä is a basin of Lake Vesijärvi (110 km 2 ) situated in southern Finland (61°10 0 N, 25°35 0 E) (Fig. 1). Enonselkä has a surface area of 26 km 2 and a mean depth 6.8 m. Originally, Enonselkä was a clear-water basin, but in the 1980s, the water quality was degraded by massive cyanobacteria blooms (Horppila et al. 1998). In the 1990s, after a largescale biomanipulation, water quality improved and harmful blooms of cyanobacteria decreased, but in the 2000s, the blooms appeared again (Keto et al. 2005). Nowadays, the concentration of TP in epilimnion is around 20 lg L -1 , but it increases above 100 lg L -1 in the hypolimnion during the summer stratification period, indicating strong internal P loading. Artificial aeration was begun in 2010 to improve the oxygen conditions in the hypolimnion and to inhibit the internal P loading. Aeration is conducted with Mixox aerators (Mixox MC-1100; Water-Eco Ltd, Kuopio, Finland) that pump surface water (3 m depth) through the thermocline (pipe diameter 110 cm) to the hypolimnion (outlet 10 m from the bottom), without breaking the thermal stratification (Lappalainen 1994). The mass flow of one Mixoxunit is 1.007 m 3 s -1 (87,000 m 3 day -1 ), and during the open water season, the pump is usually working from the onset of thermal stratification in June until the autumn circulation in August-September. In total, there are nine pumps installed in the different deeps of the basin. In the present study, the aerator situated at the Lankiluoto deep (depth 30 m, Fig. 1) was used.

Methods
The effects of aeration on temperature distribution were studied by a two-step approach. First, 3-D grids were build up based on lake topography and shape of the pump pipe, with denser grid cells around the pump (Fig. 2). Second, a 3-D model was developed for calculating the vertical temperature distribution after aeration. This numerical model was based on transient Navier-Stokes equations, continuity equations, and energy equation, which can predict 3-D transient velocity, relative pressure, temperature, and density distribution in the lake. K-Epsilon model was used to predict turbulence properties with consideration of wind-wave convective heat transfer and evaporation. The model structure and numerical simulations

Full-ported grids
The topology of the Enonselkä basin is quite complex, and it was supported by isobathic numbers. To build up 3-D grids, these measured data should be a queue on the basin surface. Depth was set as vertical (Cartesian coordinate z), while water surface was set as horizon (Cartesian coordinate x, y). MS Excel was used for horizontal 2-D lake area plots. One of the points (horizontal coordinate) on the boundary was chosen for an initial point. A FORTRAN program was developed for searching the next points. The standard of this program was the minimum distance. If the horizontal distance was the same with the other points, depth was used for comparisons. Finally, 3-D geometric grids were established by Gridgen 15.16 (licensed to Hangzhou Dianzi University). There were 34 layers of vertical grids, while horizontal grids were not restricted, and there were 7.8 9 10 6 grid cells in all.
Jacobian function in Gridgen 15.16 was carried out with all blocks (115). The value of 'y ? ' was calculated under the smallest speed of water flow near the sediment under the following equation (Wang 2014): where m is fluid viscosity (m 2 s -1 ), s x is wall shear stress (Pa), Dy is the distance between the fluid and wall, and q is water density (kg m -3 ). y ? values of the bottom and surface were calculated and they varied between 25 and 170, the average value being 89 (less than 200).

Numerical simulation model
For freshwater lakes showing salinity between 0 and 0.6, temperature 0-30°C, and pressure 0-1.8 9 10 7 Pa, the density of the water can be calculated following Chen and Millero (1986): where t is the temperature in°C. The Navier-Stokes equation is suitable for modelling water movements and heat exchanges (Pan and Eranti 2009). The velocity components in Cartesian co-ordinates for x, y, and z are named u, v, and w, respectively.
The continuity equation was expressed as follows: The momentum equations were: In addition, the heat transfer equation Q: where x, y, and z are Cartesian co-ordinates; q is density (kg m -3 ); and u, v, and w are velocity components. p is pressure (Pa); s is shear stress (Pa); and g is the acceleration of gravity; C p is the specific heat (J K -1 ), T temperature (°C), l t the turbulent viscosity (Pa s), and _ q the generation of thermal energy per unit volume (W m -3 ).U is the viscous dissipation function. C is expressed as where k is heat conductivity and Pr is the turbulent Prandtl number.
At the surface boundary, wind stress can be treated as a source, which affects horizontal movements and disturbs stratified water layers (Laval and Imberger 2003). Wind stress on the water surface can be calculated by: where q a is the density of air (1.2 kg m -3 ), C D is coefficient of drag, and u w is the wind velocity (m s -1 ). The first iteration of numerical diffusion in simulations was controlled by wind stress, and transient (every 10 min) wind data were used. In accordance with velocity, wind direction needs to be decomposed into x and y directions, which can be expressed by the following: Water flow near the surface was calculated by a K-x model. Wind-wave induced turbulence (core areas) as well as the turbulence induced by the pump was evaluated by a K-epsilon model.
The hybrid scheme (Versteeg and Malalasekera 2010) was used for spatial discretization of momentum equations, heat transfer equation, as well as for K-epsilon equations. Full implicit scheme was used for temporal discretization. The continuity equation was not solved alone, but was coupled with momentum equations using the SIMPLE scheme to get pressure (Pan 2000;Smith et al. 2000).
Water temperature is affected by solar energy and thermal radiation from the water to the air. Water can absorb heat from global irradiance when the value is [120 W m -2 (Badescu 2016). The average of global irradiance was considered as a heat generator tiled on the lake surface. The net emission of thermal radiation from water surface (cal cm -2 day -1 ), was calculated as follows: 11 9 water temperature (K)-air temperature (K) (Hutchinson 1957;Wetzel 2001).
These discrete equations and solar energy equations were solved by a FORTRAN program. Every iterated step took 75 ms with Core-i5-6500 and 12 central processing unit (CPU). Therefore, 1 day simulation took 108 min., and 5.4 and 27 h were used for modelling 3 and 15 d aeration periods, respectively.

Boundary conditions and collection of field data
To study the aeration effects, the vertical temperature distribution of Enonselkä on August 4th 2009 without artificial aeration was set as an original water temperature. Data on wind velocity and direction (Fig. 3) were obtained from the Laune meteorological station (Finnish Meteorological Institute), which is situated 3 km from the Enonselkä basin (measurements every 10 min).
For global irradiation, values[120 W m -2 were used as a heat producer. The value was 257.407 W m -2 during 30 July -1 August and 307.386 W m -2 during August 7-22, respectively ( Table 1). The data for the solar energy were obtained from the Laune meteorological station. The inflow of rivers to Enonselkä is small; the mean inflow of the four largest brooks together is only 1.1 m 3 s -1 and their effect was thus excluded ( Table 1). The relative pressure of Enonselkä basin outflow is 0 Pa, while the mass flow is free, because the water in the Enonselkä basin is in connection with the other parts of the lake (Fig. 1).
The pump was set as a cylinder tube; volume flow of this pump was 1.007 m 3 s -1 (87,000 m 3 day -1 ). The onsets of aeration were set for 7:00 a.m. 30 July and 0:00 7 August, respectively, and these times were also when the model started. The lake surface was set as free face, while the lake bottom, coastline, and island boundaries were set as no slip and heat insulation wall. The reference pressure of the lake was 1 atm. Buoyancy forces and eccentric forces generated by Earth's rotation were considered.
The vertical temperature distribution of Enonselkä on August 4th 2009 (Fig. 4) was measured with a YSI-6600V2 sonde (YSI Corporation, Yellow Springs, OH, USA) and was set as a boundary condition. It presented conditions during summer stratification. This is the period when aeration is especially used due to the deterioration of oxygen conditions in the hypolimnion during the summer. Water temperature in the epilimnion was 18-19°C, thermocline was situated at 6-8 m depth, and in the hypolimnion, temperature decreased from 15°C at 10 m to 9°C at 30 m depth (Fig. 4).
The effect of aeration on temperature distribution and model performance was studied using the model in predicting the temperature distribution of Enonselkä on two dates, 1 August and 22 August 2013. For research purposes, the aerator was not in use during 14-28 July and was started in 29 July. Thus, on 1 August, the aerator had been functioning for 3 days. Another break in aeration took place from 2 August to 6 August. On 22 August, the pump had thus been functioning for 15 days. The instantaneous wind data from July 30 to August 1, 2013 and from 7 August to 22 August 2013 were used as an input. From the lake, the observed temperature profiles on 1 and 22 August were measured with the YSI-6600V2 sonde.
The relationship of the predicted temperatures with the observed ones was studied with linear regression analysis. In addition, the residuals of the model were plotted against the observed values, and the relationship between residuals and observed values was analysed with linear regression. The analyses were conducted separately for epilimnion and hypolimnion, because only hypolimnion was assumed to be affected by the aeration. Similarly, the relationship of model residuals with water temperature and depth in the epilimnion and hypolimnion was studied with regression analysis. The data sets were ln-transformed before the analyses.

Results
On 1 August 2013, after 3 days of aeration, the temperature stratification was steeper than in the boundary conditions (Fig. 4). Thermocline was situated at 9 m depth and temperature in the hypolimnion was lower than in boundary conditions. The modelled temperature in the surface (18.7°C) was 0.8°C lower than the observed temperature (19.5°C) (Fig. 4). At 5 m depth, the model gave a prediction identical to the measured temperature (difference \0.1°C), while at 9 m depth, the model underestimated water temperature by 1.7°C. The difference between the modelled and observed temperatures was largest at the lower part of the thermocline (10 m depth) where it was 3.6°C (modelled 12.6°C, observed 16.2°C). At the 15-20 m depth in the hypolimnion, the modelled and observed temperatures were close to each other, whereas below 20 m, the model overestimated water temperature by 0.2-0.5°C (Fig. 4). The observed and modelled temperatures had a tight relationship (F 1,20 = 532.65, r 2 = 0.964, p \ 0.0001) (Fig. 5).
On 22 August, after 15 days of aeration, the observed temperature profile was closer to boundary conditions (Fig. 3). Temperature in the hypolimnion was higher than at the boundary conditions. The modelled temperature in the metalimnion was 0.4-0.9°C higher than the observed temperature (Fig. 4). In the thermocline, the modelled and observed values were close to each other with the exception that the modelled thermocline extended deeper (14 m) than the observed one (12 m). In the hypolimnion, the model underestimated the observed temperature (11-13°C), the model error values decreasing towards the bottom. In all, however, model prediction closely followed the observed temperature (F 1,19 = 209.30, r 2 = 0.921, p \ 0.0001) (Fig. 5). Model residual did not depend on water temperature in the epilimnion (F 1,5 = 3.29, r 2 = 0.397, p = 0.130), or in the hypolimnion (F 1,11 = 4.60, r 2 = 0.295, p = 0.06) (Fig. 6). The effect of water depth on model residual was similar to 1 August; residual increased together with depth in the epilimnion, but the regression was not significant (F 1,5 = 6.36, r 2 = 0.560, p = 0.053) (Fig. 7). In the hypolimnion, residual decreased significantly with depth (F 1,9 = 532.89, r 2 = 0.987, p \ 0.0001) (Fig. 7).

Discussion
The model predictions were close to the observed temperatures. On both study dates, the average difference between the modelled and observed temperature was 0.6°C. The average percentage model error (all depths pooled) was 4.0% on 1 August and 3.7% on 22 August. The aerator had a clear effect on water temperature in the hypolimnion. On 1 August, when the aeration had been functioning for 3 days, temperature in the hypolimnion was 11.1-12.5°C. On 22 August, when the aerator had been functioning for 15 days, temperature in the hypolimnion was 14.8-15.3°C. The change was not caused by weather conditions, because at the same time, the average temperature in the epilimnion decreased by 0.3°C. Thus, the increased temperature in the hypolimnion was due to the warm epilimnetic water transferred downward by the aerator. Such temperature effect of hypolimnetic aeration has previously been described in numerous studies (Soltero et al. 1994;Beutel et al. 2007;Gantzer et al. 2009).
In the epilimnion, the accuracy of the model prediction was dependent on the difference between prevailing temperature and boundary conditions. Temperature in the epilimnion was slightly underestimated on 1 August when observed temperature was higher than in boundary conditions and overestimated on 22 August when observed temperature was lower than in boundary conditions. However, because the aim of the model was especially to study the effects of aeration in the hypolimnion, and mass flow between the hypolimnion and epilimnion is weak, the effect of the boundary conditions did not affect the feasibility of the model.
In the epilimnion, model residual was not dependent on temperature or depth, but in the hypolimnion, the residual decreased with decreasing temperature (1 August) and decreased with increasing depth (both dates). These two phenomena are connected as temperature decreases with depth, but the effect of depth was stronger than the effect of temperature. The model was thus not able to precisely describe the vertical distribution of heat in the hypolimnion. On 1 August, the model predicted a homogenous temperature profile in the hypolimnion, while the observed temperature decreased moderately from the thermocline to the bottom. This was probably because the effect of sediment was not included in the model. In summer, the sediment heat flux is from the sediment to the water for the mixed layer and from the water to the sediment below the mixed layer (Fang and Stefan 1996). Thus, in summer in the stratifying areas, the sediment reduces the temperature of the overlying water (Lindenschmidt and Hamblin 1997), which resulted in the slight overestimation of the nearbottom water temperature by the model. According to Fang and Stefan (1996), in lake with the morphometry of Enonselkä, the sediment effect on bottom water temperature is 1-2°C, which was in accordance with the present study. On 22 August, the modelled and observed temperatures near the bottom were identical, demonstrating that the heat transfer by the aerator had masked the effect of sediment on temperature and suggested that exclusion of sediment heat from the model does not cause considerable error unless very short-term effects of aeration are studied. In all, the model successfully described the effects of the aerator on the lake's temperature profile. The results confirm the validity of the applied transient Navier-Stokes equations and demonstrated that iteration every 0.5 s allowed a timely calculation of all the instantaneous parameters.
creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.