Ensemble Modeling of the Impact of Climate Warming and Increased Frequency of Extreme Climatic Events on the Thermal Characteristics of a Sub-Tropical Lake

: Lake ecosystems are impacted by changes in climatic conditions. Climate changes forecasted to occur are reﬂected in models by slow gradual changes over extended periods of time. Output from weather generators, on the other hand, can simulate short-term extreme conditions and weather patterns. In order to evaluate the likely impact of climate changes on a large sub-tropical lake, speciﬁcally the thermal regime of the lake, we constructed climate scenarios using a weather generator. The 30-year scenarios included no change in climate conditions, a gradual change, increased frequency of heat waves and a merging of the latter two. The projected impact on the lake’s physical properties was evaluated using an ensemble of 1-D hydrodynamic lake models. The gradual increase scenario had the largest impact on annual temperatures and stratiﬁcation period; however, increased heat waves had a large e ﬀ ect on the summer lake conditions and introduced a larger degree of variability in water temperature. The use of the ensemble of models resulted in variability in the projected impacts; yet, the large degree of similarity between projected trends and patterns increased conﬁdence in the results. The projected e ﬀ ect the heat waves will have on the lake conditions highlights the need to include heat waves in climate studies and the need for impact studies in order to better understand possible consequences for lake ecosystems.


Introduction
Climate models are widely used by scientists to explore the present, as well as climate change scenarios, and serve as a basis for policy debates and decisions [1]. In an attempt to apply global climate change models (GCMs) to spatially limited areas and use higher resolution models, a number of approaches have developed, which include regional climate models (RCMs), statistical downscaling and weather generators [2,3]. Each of these approaches provides different benefits, but they suffer from various limitations, and their efficiency depends on the objectives of the study at hand.
One of the shortcomings of climate models, regardless of the spatial and temporal extent and resolution, is the monotonous change over time in the model variables. As a result, changes forecasted to occur over time are represented by slow gradual changes over extended periods of time. These An increase in the frequency and duration of heat spells • An increase in the summer heat load, especially during August, the hottest month.
Lake Kinneret, the largest natural freshwater lake in Israel and the Middle East, is a key source of drinking water, located at 210 m below sea level (−210 asl), and is in a region that experiences high summer temperatures. It has been reported that the number of days with temperatures above 39 • C has increased by 5.6 days per decade and the upper 95th percentile of the daily maximum temperature for the period 1975-2010 was 37.4 • C [25].
The increased frequency of extreme events and overall warming trend is likely to have multiple effects on a lake ecosystem such as Lake Kinneret. The initial effect will be on the physical characteristics as they are impacted directly by any change in the heat fluxes. In fact, the annual maximum water temperature (at Z = 1 m) of Lake Kinneret has already demonstrated a long-term (1970-2019) trend of increasing surface water temperature at a rate of 0.4 • C per decade ( Figure 1). This rate of increase is higher than previously reported for the epilimnion over the period 1969-2008 [11]. In addition to the long-term warming, the lake has been exposed to increased frequencies of extreme heatwaves (three-day episodes with a maximum temperature above 35 • C), the most prominent to date occurring in 2010. Based on the long-term weekly profiles conducted since January, 1969, at a station located at the deepest section of the lake (Station A, Figure 2), the highest water temperature recorded (31.2 • C) occurred on 22 August, 2010. Data collected from a higher resolution thermistor chain indicated that the maximum water temperature was 34.32 • C (recorded on 20 August, 2010) and that 68% of the water temperatures above 32 • C recorded since January, 2000, occurred during the heat wave of August, 2010. Short-or long-term changes in the heat budget and thermal characteristics can percolate through the ecosystem by altering process rates. For example, key metabolic rates of all biological organisms are temperature dependent. Thus, increased water temperatures will lead to an increase in metabolic rates and energy demand [26]. As a consequence, the observed climatic trends are likely to have a notable impact on the lake ecosystem.  Figure 2). The cumulative distribution (B) includes all available measurements collected between 2010-2018 by an automatic profiler, mounted with a Manta 35+ multiprobe, 4-6 times daily. Profiler data were not available for most of August, 2019, therefore 2019 is not included. The profiler is mounted on a fixed station located at Stn. A.
In order to evaluate the likely impact of climate change on the Lake Kinneret ecosystem we constructed climate scenarios using a vector-autoregressive weather generator [27]. The climate scenarios served as input into three 1-dimensional lake hydrodynamical models, which allowed us to conduct ensemble simulations of the impact on the lake thermal conditions. Based on the ensemble of models, we evaluated the projected impact of climate warming and increased occurrence of heat waves on the lake. The current study is unique in the merging of lake-ensemble modeling with lake-specific climate change scenarios accounting for not only a general trend of increasing temperature but also the increased occurrence of extreme heat events.

Study Site
Lake Kinneret is a sub-tropical monomictic lake located in the northern part of Israel ( Figure 2). The 45 m deep lake (at maximum lake level) stratifies during the spring and remains stratified through the summer and fall to complete mixing typically in December or January. Two key climate characteristics impact the lake's heat budget: high air temperatures throughout the year and a strong warm afternoon breeze during summer months [28]. Based on long-term data measured at a meteorological station >3 km south of the lake, the annual average air temperature is relatively high (21 • C) with maximum summer air temperatures exceeding 36 • C [28]. Seasonally, the summers (May-September) are dry, and the winters are short with a mean rainfall of 380 mm and maximum daily air temperatures ranging typically between 15-20 • C.

Weather Generator
The weather generator, a vector-autoregressive weather generator (VG), uses at its core a vector-autoregressive (VAR) process [27]. A VAR process separates a multivariate time series into deterministic and random parts. The variables at each day are assumed to be linearly dependent on the values of all variables from the previous several days. Estimating the parameters of the deterministic part and drawing random numbers that correspond to the subsequent residual part, one can generate arbitrarily long synthetic time series. As such, time series exhibit normal marginals, and transformations based on quantile-mappings are used to make them resemble the observed reality ( Figure 3). The data used as the basis for VG data generation were collected at an on-lake site located in the north-west part of the lake, approximately 1.5 km from shore (depending on lake level). The meteorological station, located at this site since 1996, recorded air temperature, relative humidity, global short-wave radiation, incident long-wave radiation and wind speed and direction [28]. The time span used for this study includes 1997 to 2007 (11 full years).
For generating scenarios, disturbances are added during time series generation. These disturbances were applied to all variables, in accordance to their correlation structure. This means that a change in, e.g., air temperature is propagated to other variables that are seen to reflect changes in air temperature. This design was chosen deliberately so that (1) scenarios can be defined easily using just one variable, and (2) those scenarios are plausible in terms of their statistical dependencies.
Compared to [27], in which the weather generator was calibrated on a data set measured in Central Europe, a number of changes were made in order to represent the different climate. Air temperatures at Lake Kinneret are not symmetric to a seasonal mean, but skewed towards higher values. As these higher values are related to heat waves, we decided to use the Gumbel instead of the normal distribution. Heat waves at Lake Kinneret are often related to an obstruction of the Mediterranean Sea Breeze (MSB) by a stable inversion layer [29]. These heat waves do not manifest as solitary increases of air temperature; they are also accompanied by a decrease in humidity and wind speed. Reproducing these multivariate patterns was achieved by VG through generating random sequences of low relative humidity (25%) events during the months June through September. VG was updated to accept changes to more than one variable in order to produce combined scenarios that contain modifications to both air temperature and relative humidity. For each changed variable, the influence on the other variables was calculated and added separately as described for single variable changes in [27]. Thus, gradual change and heat waves can be seen as different drivers, each having an influence on the full set of simulated variables.

Hydrodynamic Models
The climate scenarios were applied to a three-model ensemble that included the following 1-D hydrodynamic models: the dynamics reservoir simulation model (DYRESM), the general lake model (GLM) and the general oceanographic turbulence model (GOTM). DYRESM has been used extensively to simulate the thermal dynamics of lakes [30]. It has been calibrated and tested for Lake Kinneret [31,32], used to examine changes to the lake thermal structure as a result of lake level changes [11], and has been used as a hydrodynamic driver for the application of the biochemical Computational Aquatic Ecosystem Dynamics Model (CAEDYM) to Lake Kinneret [33]. GLM is a more recently developed model similar in process description to DYRESM [34]. It has been applied to date mainly to temperate lakes [35,36] though it was used for a large scale comparison between lakes, including Lake Kinneret [37]. GLM combines fluxes of mass and energy with a Lagrangian layer structure that adapts to changes in vertical gradients. Detailed GLM formulations of the numerical routines can be found in [34]. For the application of GOTM we used the lake-adapted version of the model. While originally developed for marine applications [38] it has since been applied to lake systems [39][40][41]. Details about the model can be found in Burchard et al. [42].
As this is the first application of GOTM and a calibrated version of GLM to Lake Kinneret, we provide a description of the calibration and validation of both models. We do not provide results of the DYRESM calibration and validation to Lake Kinneret in the current paper, but rather refer the reader to the study by Gal and colleagues [32] and their testing of the model on Lake Kinneret.

GLM Calibration and Validation
We calibrated and validated GLM to Lake Kinneret based on data collected as part of the ongoing monitoring of Lake Kinneret [43]. The calibration-validation period extended between 1 January, 1997, and 31 December, 2003, where the data from January 1997 to December 1999 was used for calibration and the remaining data for validation. The model was set up to allow a layer thickness of 0.24-1.54 m and an hourly time step. The temperature data were collected at a weekly resolution at 1 m intervals from the surface to a depth of 40 m, thereby providing information on the thermal structure of the lake. All data were collected at Station A (Stn. A). The inflows input data consisted of discharges (volume per day), temperature and salinity, whereas the outflows included only discharge. There were a total of six inflows of which approximately 70% of the incoming water was accounted for by the Jordan River [32]. Inflows and withdrawals for all three models were based on water balances and annual reports published by Mekorot, Israel's national water company.
Additional inputs included meteorological variables (short wave radiation, long wave radiation, air temperature, wind speed and relative humidity), lake bathymetry and initial profiles of water temperature and salinity and lake level at day 1 of the simulation. Time varying input data are provided at daily resolution with the exception of the meteorological data input, which are provided at hourly resolution. We refer the reader to [32] for a more detailed description of the model input data.
The model was calibrated based on a three-year period using the genetic algorithm (GA) approach (see below) and then tested, without performing changes to the parameters, on an additional four years. The model time step was hourly but the output was daily. The output was compared to temperature data for the same dates and depths. Model performance was assessed by comparing the surface (0-10 m) and bottom (30-40 m) layer temperatures using the Nash-Sutcliffe efficiency index (NSE), R 2 and root mean square error (RMSE) statistics.

Genetic Algorithm
Genetic algorithms (GAs) use an initial population of random possible solutions and determine which of these possible solutions are the strongest using a fitness function. The fitness function was defined, in this case, as the minimum value of sum of square errors (SSE) between predicted values and observed values. In this study, the fitness function (Equation (1)) was based on the temperature from the bottom layer (30-40 m) and from the surface layer (0-10 m) during the spring and fall when stratification is developing or eroding, as these are the more challenging periods in terms of simulating water temperature. The individuals with the strongest fitness are then "bred" and produce "offspring". After a crossover is performed, mutation takes place by randomly changing the new offspring. The offspring are again evaluated to distinguish strong from weak. In each successive generation, new parents are selected and offspring are produced. This process continues until the process converges on an acceptably strong solution or reaches one of the stopping criteria. We defined the fitness function as: Fitness function = min (sum(e1+e2)) (1) where e1=1-NS_bottom, e2=1-NS_strat and NS is the Nash-Sutcliffe coefficient value (NSE); bottom represents the 30-40 m layer of the water column and strat represents the surface layer (0-10 m) temperature only during the periods of the spring development of the stratification (April-June) and the fall erosion of the stratification (October-December). The min refers to the minimum value of the vector of values calculated for the entire range of solutions. The following limitations were imposed on the fitness functions in order to avoid negative values: If NS_strat<0 then NS_strat=0 If NS_bottom<−1, e1=|NS_bottom| The population size at each generation was set at 100 individuals, who were presented as real numbers. Stopping criteria were defined as the conditions in which the cumulative change in the fitness function value is less than the function tolerance (10 −6 ) or if the total allocated number of generations is reached. Solutions generated using GA generally show a rapid evolution and an immediate increase in fitness scores in early generations, followed by very small changes in later generations; therefore, it seems inefficient to set the generation number to an arbitrarily high value. The optimal solutions were searched using the Matlab genetic algorithm optimization function.
The parameters that have the greatest impact on the evolutionary process and the rate of convergence toward the optimal solution include proportion of crossover, proportion of mutation and the presence of elites, the best individuals that are shifted to the next generation unchanged.
Crossover improves the genetic selection and combination, whereas mutations are generally damaging rather than beneficial to the optimization process, as they are in nature. However, mutations provide a mechanism for the search to jump out of a local optimum point. The proportion of crossover, mutation and elites in the population are normally determined by trial and error, as there is no universal method to determine the optimal GA operators' values. We used the Matlab GA function default values of 80% crossover, 10% elite and 10% mutation, which we have found in the past to provide satisfactory results [44].
Since the GA does not necessarily find the global optimum of the fitness function, but may converge to a local optimum, it must operate for a number of runs and the best solution will be determined according to the minimum values of the fitness function [44]. Therefore, in this study the calibration procedure was repeated for 10 GA runs to achieve the best solution. Preliminary trials indicated that having a larger number of runs does not improve the results of the GA calibration.
The GA-based calibration was applied to the parameters in the GLM parameter file. A total of 13 parameters were included in the calibration process with initial values based on the software default values ( Table 1). The ranges of acceptable values for each parameter were in most cases 20% above and below the default value, though for some parameters we extended the range based on preliminary testing we conducted.

GOTM Calibration and Validation
The model was set to simulate the period January, 2010, to September, 2019, where the first three years were defined as a spin-up period, which is often used in order to allow the model to reach stability prior to evaluating performance. The calibration was conducted only on the remaining period. Input data included high frequency meteorological variables (air temperature, relative humidity, air pressure, wind speed, shortwave radiation) collected at Stn. A. The simulation's time steps were set to an hourly frequency with a vertical resolution of 50 equal-depth layers.
Model calibration was performed using an automatic calibration program-ACPy [45]-to adjust a set of parameter values to provide the most accurate results. A total of seven parameters were included in the calibration process ( Table 2). The ACPy program used for calibration is based on a differential evolution algorithm that compares simulated and measured water temperature and calculates the log likelihood function. ACPy is based on an iterative process, where each iteration runs the model with different values for the pre-selected calibrated parameters and compares the results to observed data using the log likelihood function [46]. Before initiating the process, the user defines the parameters to be calibrated and their ranges of values. During the iterations, the parameter values converge towards values with a greater likelihood, until achieving a set of values with the best likelihood. ACPy was set to conduct a total of 5000 iterations of the calibration process. For the calibration, we used water temperature profiles measured by an auto profiler device with a Manta +35 multiprobe (Eureka Water Probes, Austin Tx, USA) located at Stn. A. that profiles the water column four to six times daily.

Ensemble Scenario Testing
We conducted a series of scenario simulations in order to determine whether the projected climate change and observed changes (i.e., [25]) in meteorological characteristics will have an impact on the thermal characteristics of Lake Kinneret. The simulations were composed of climate scenarios, generated by VG, based on predicted and observed long-term changes using an ensemble of DYRESM, GLM and GOTM. We examined four different climate change scenarios varying in the characteristics of the meteorological input data. The four scenarios included (Table 3): 1.
Unchanged-current meteorological characteristics remained unchanged. This scenario served as a reference for the other scenarios.

2.
Linear Increase (LI)-a gradual linear increase in air temperature over the entire period of the simulation. Temperature increased by a total of 2 • C over the period of the simulation.

3.
Heat waves (HW)-increased frequency of extreme heat events during the summer similar to the changes described by Ziv and colleagues [25]. 4.
HW-LI-Merging of increased frequency of extreme events (scenario 3) and a gradual increase (scenario 2). This scenario is considered the most realistic.
The sequence of random numbers generated within VG depends only on the realization's index, making each realization comparable throughout the scenarios (e.g., in the nth realization, heat waves occur with the same timing and length in both the HW and HW-LI scenarios). For each scenario we conducted 1000 30-year realizations (thus, in total 4000 simulations for each model) during which there was no interannual variability in the input data, with the exception of the meteorological data. Inflow and outflow/withdrawal data were defined as for the calibration process for all the models. The inflows input data consisted of volume, temperature and salinity, whereas the outflows included only volumes, both at a daily resolution. In both cases inputs were based on actual data from the year 2002, which was considered a typical year.
All three models were configured identically for the scenario testing. Simulations were initiated on 1 January, 2010, with an hourly time step and daily output. The same input data were used for all models accounting only for differences in input file structures. Parameter values were based on the calibrated values. The extinction coefficient value (K d ) was set at 0.6 m −1 in both DYRESM and GLM.
The results we present are based on all 1000 realizations of each of the four scenarios. We evaluated the impact of the various scenarios by analyzing a number of metrics focusing mainly on the later period of the simulations in order to examine the long-term effects. The initial five years of the simulations were not used for any part of the analyses, in order to eliminate the effect of initial conditions. Temporally, we conducted long-term annual and seasonal analysis and shorter-term analyses of the final years of the simulations. The main focus of the analysis was on the upper layer of the water column (0-10 m), which is where a majority of the biological activity occurs and is the main portion of the lake water column most expected to be impacted by climate change. Metrics included in the analyses were annual mean temperatures and divergence from the base scenario (termed "Unchanged"), summer water temperatures, as well as duration of stratification as increased heating is likely to increase stratification stability. The stratification period was defined by the number of days in a year in which the difference between the mean 0-10 m temperature and the mean 30-40 m temperature was greater than 2 • C.

Calibration Results
Results of the GLM calibration for the lake indicated the model was able to accurately simulate the thermal characteristics of the surface layer of the lake with a Nash-Sutcliffe efficiency (NSE) of 0.98 over the entire simulation period (Figure 4, Table 4). During the period of stratification set-up and break down, the NSE was still high, with a value of 0.86. The model, as expected, did not perform as well in simulating the dynamics of the bottom layer water temperature, resulting in an NSE of −3.30. The lack of accuracy was mainly due to the reduced range of simulated temperatures (13.52-16.96 • C) in relation to a somewhat wider range of measured values (13.74-17.28 • C). The GOTM calibration provided a more accurate representation of the lake data, with NSE values of 0.98 and 0.18 for the surface and bottom layers, respectively, and low RMSE values for both layers (Table 4). Although both models underestimated water temperature in the bottom layer with a mean difference of 1.07 • C for GLM and −0.15 • C for GOTM, the timing and amplitude of the surface layer seasonal heating were accurate, as was the annual water column mixing period and the seasonal and interannual trends in the bottom layer (Figures 4 and 5).

Scenario Results
Analysis of the changes in temperature over the 30-year period of the scenarios resulted in minor changes between the altered climate scenarios and the Unchanged scenario. The mean increase in temperature of the upper layer in the heat wave (HW) scenario, in relation to the Unchanged scenario, ranged between 0.00-0.12 • C over the ensemble of models, though there was, as expected, a large degree of variation in the results between the realizations and over the ensemble of models ( Figure 6). As the average was calculated over the entire year, this outcome is not surprising. The increase in the Linear Increase scenario ranged between 0.53-0.55 • C, whereas the combined HW-LI scenario resulted in a mean increase of 0.54-0.63 • C. The overall range of changes in the two extreme scenarios represented an increase of only 1.3%-2.7% (Figure 7). The increase in the temperature of the bottom layer (30-40 m) was moderate and ranged between 0.21-0.51 • C in the two extreme scenarios, representing a 2.5%-8.4% increase in water temperature in relation to the Unchanged scenario. In addition to the mild mean increase in water temperature over the entire period, changes in the 90th percentile of water temperature in the upper 0-10 m of the water column occurred only in the two scenarios that included a gradual increase in temperature ( Figure 8B,C). Although there was no difference between the HW and Unchanged scenarios in the value of the 90th percentile at the end of the simulation, the differences over the entire simulations (and ensemble of models) ranged from −0.19 • C to 0.84 • C. The differences between the LI and Unchanged scenarios were larger and ranged between −0.01 • C and 1.12 • C, reaching values of 0.68 • C to 0.81 • C at the end of the simulations. The HW-LI scenarios resulted in a larger range of increase in the 90th percentiles in relation to the Unchanged scenario, with differences ranging between −0.09 • C and 1.67 • C, reaching a difference of 0.81 • C during the last month of the GOTM simulation. The addition of the HW conditions to the LI scenarios introduced substantial variation in the range of values (Figure 8), highlighting the impact of heat waves on the extreme summer water temperatures. The impact of the increased frequency of heat waves in the HW scenarios on the water temperature is particularly noticeable during the summer months. As a result, the mean values of the HW scenarios are higher than the median values, reaching a difference of 0.1 • C in the GLM simulations ( Figure 9). Results of the HW-LI scenario display an extended distribution of temperatures, reaching temperatures up to 1.6 • C higher than the LI scenario (based on averages of all realizations). Though the probability of these extreme temperatures is low, the pattern is nevertheless consistent over all three models and similar to the pattern observed during August 2010. It is interesting to note that the results of the final years of the LI scenario have lower mean and median temperatures than those measured in August 2010, though the extreme temperatures are similar. The 2010 mean and median temperatures are, however, notably higher than those measured in recent decades (Figure 9). Occurrence of temperatures at least as extreme as 2010 is most probable in the HW-LI scenario, thus highlighting the influence of the heat waves. It is noteworthy that the Unchanged results are similar or cooler than the lake data since 2000. This is because VG was calibrated on the cooler 1997-2007 period.    The LI scenario resulted in a gradual increase in water temperature; however, the impact of the heat waves introduced by the HW scenarios resulted in short periods of high temperatures. This was particularly evident during the early summer months (Figure 10). The later part of the summer was impacted to a larger extent by the gradual increase in temperatures, and in both cases resulted in a warmer, more stable epilimnion. This was clearly evident in the combined HW-LI scenario. Heat waves combined with a gradual increase are followed by a longer-lasting increase in summer epilimnetic temperatures than under present-day temperatures. A consequence of the stabilizing of the epilimnion was the increase in the stratified period ( Figure 11). For all three models there was a significant difference (Kruskal-Wallis one-way ANOVA on Ranks and Tukey, p < 0.05) in the length of the stratified period between all scenarios, with the exceptions of HW vs. Unchanged and HW-LI vs. LI. It is noteworthy that the observed long-term mean length of the stratification period based on data collected between 1969 and 2008 was approximately 286 days and remained stable over this period [11]. This value is comparable to the results of the GLM and DYRESM Unchanged scenario, serving as an indirect validation of the model. The simulated GOTM stratification period in the Unchanged scenario is shorter than the period simulated by GLM and DYRESM.

Discussion
The main goal of this study was to evaluate the likely impact of the occurring and projected climate changes in the Eastern-Mediterranean, and in particularly Northern Israel, on a sub-tropical deep lake based on an ensemble of 1-D hydrodynamic lake models. It was not our intention to provide an extensive comparison between the models and a review of the possible differences, but rather to take advantage of the variations between them. This study represents one of only a limited number of lake ensemble modeling studies or variation of ensemble modeling approaches [12,47,48]. The use of an ensemble and the range of differences between the models offer a wider view of the possible impact of the increase in frequency and intensity of heat waves occurring in this region. We further introduced variability between the models by applying three models that had been calibrated based on different periods. Thus, given the large potential source of uncertainty between the models, the observed similarity between them in response to the climate forcing scenarios strengthens our conclusions.
The projected increase in air temperature (2 • C over a period of 30 years), as represented in our LI scenario, had a clear long-term effect on lake water temperature. The ensemble results of the LI scenario indicated an annual mean warming of the lake, an increase in the maximum values, warming of the summer water temperatures and an extended stratification period. In contrast, the HW scenario, in which we introduced increased frequency of extreme heat events during the summer similar to the changes described by [25], resulted in a number of effects. These included a moderate annual increase in water temperatures but mainly increased daily and summer water temperatures as well as increasing the variability of the maximum water temperatures. The projected extension in the length of the stratified period, though there was only a moderate increase in water temperature, is in line with global long-term data and additional modeling studies [49]. A daily increase in water temperature occurred during extreme heating events, but when merged with the LI scenarios the extreme events had a prolonged impact that extended through to the end of the summer. In this detail, the influence of a combination of increased climatic mean and variability is more than an addition of its individual influences. This exemplifies the merit of our four-scenario approach, and points to possible threats, given the trends described by Ziv and colleagues [25]. On an annual scale, the merging of the extreme events with the gradual increase in air temperature had, however, only a minor impact. This was a consequence of the definition of the HW scenarios, in which extreme heat events were limited to the summer period only.
The overall annual increase in water temperatures and large increase in mean temperature of the upper layer could have major implications on the ecosystem. Stratified lakes are particularly sensitive to changes in the weather as the surface layer integrates the effects of changes in the surface heat fluxes driven by solar radiation, air temperature and wind speed, and retains the heat until water column mixing. However, in the Mediterranean and other sub-tropical regions, the warmer temperatures may also promote different effects than those expected in temperate regions, because they have higher annual water temperatures and higher seasonal fluctuations in water level, but smaller seasonal changes in light and temperature, in comparison to temperate lakes. It has been shown, based on long-term temperature data from 26 lakes, that deep and warm lakes have experienced the largest changes in lake stratification as a result of climate warming, and that the magnitude of change in lake stratification is primarily controlled by lake morphometry (mean depth, surface area and volume) [49]. The potential effects associated with warming and changes in stratification would be further aggravated by the occurrence of the heat waves and especially the increasing extreme temperatures. An increase in water temperatures would increase metabolic rates such as respiration and decomposition rates [26,50] and possibly lead to changes in the zooplankton community [51,52]. Although the impact of warming on phytoplankton can be dependent on trophic interactions [53] and the likely outcome for Lake Kinneret is not clear, there is strong evidence linking climate change and warming to the proliferation of cyanobacteria and intensification of harmful algal blooms [7,[54][55][56]. Furthermore, long-term global data show that the impact of lake warming on metabolic processes is likely to be the strongest at low latitudes and low elevation lakes [10]. In addition, there is concern as to the impact on water quality [57]. Thus, increased warming and frequency of extreme heat events is likely to have a major impact on the lake's ecosystem. Given the wide range of ecosystem services provided by the lake, and its unique role as a major water and recreation resource in Israel, the likely changes could have dramatic and disastrous consequences.
As the basis for the use of the ensemble of models, we presented the calibration results of two of the models, as this was their first implementation for Lake Kinneret. We presented a level 1 or "state" validation [58] of the models and show that they accurately capture the thermal dynamics observed in the lake including epilimnetic water temperature, seasonal stratification patterns and the length of the stratification period that matches previously reported information [11]. In both cases, however, the simulation accuracy of the surface layer is considerably higher than for the lower layer of the water column. This is similar to the results of the application of DYRESM to the lake [32] and is most likely a result of underestimated mixing down of heat, possibly due to a 1-dimensional simulation of highly spatially variable processes and wind forcing on the lake [59]. The lack of simulated accuracy of the lower layer of the water column does not, however, affect the key conclusions of this study, as the upper layer is the biologically active layer which will predominantly be affected by the changing climate conditions and increased heat waves.
Recent studies have documented the climate changes that have occurred in the Eastern-Mediterranean region and that are projected to occur into the future and their likely impact on the environment including Lake Kinneret. Reports have, in particular, highlighted the strengthening and lengthening of hot periods and changes in their seasonality [28,60,61]. Our analysis of the long-term monitoring data from the lake show that climate warming has led to a clear trend in increasing water temperature since the 1970s (Figure 9). During this period the mean and median temperatures rose from 28.01 • C and 27.95 • C during the 1970s to 29.76 • C and 29.65 • C during the second decade of the current century, with an increasing trend of 0.4 • C per decade. The long-term changes that have occurred in the thermal characteristics of the lake since the 1970s have been linked to an increase in spring and summer air temperature and lake level [11]. These observed changes in the region are having a range of effects on Lake Kinneret and its watershed, including a decline in runoff and stream flow into the lake, a decline in lake level and an increase in salinity [62][63][64].
Projections have suggested strengthening of the reported climate trends. The projections are based, however, on global- [61] or regional-scale and downscaled climate models of the region [60,65]. Although these models typically reflect the observed trends accurately, they represent relatively low spatial scale resolution and hence lack accuracy for localized trends [3]. The use of weather generators based on localized data circumvents this weakness and also allows for better representation of the local climate conditions and trends. We applied a weather generator used in the past for Central European sites and adapted it to the Lake Kinneret conditions. The use of the weather generator allowed us to introduce heat events with their observed characteristics and trends into our climate forcing data, which would not be possible with large scale climate models.
This study is one of the only studies applying ensemble modeling in an attempt to forecast the likely impact of climate change on a lake ecosystem [45,66]. It is the first, as far we are aware, to project the likely impact of climate change and increased frequency of heat waves on a sub-tropical, deep lake. Lake Kinneret, like many deep sub-tropical lakes, is monomictic, with an extended stratified season and an anoxic and nutrient rich hypolimnion, with a resetting of the ecosystem occurring with the annual mixing of the water column. The mixing is a key process for the physics, chemistry and ecology of the lake. Changes to the mixing and stratification sequence as a result of climate warming and extreme heating events could have dire consequences on the ecosystem and the services it provides.