Simulations of future changes in thermal structure of Lake Erken: proof of concept for ISIMIP2b lake sector local simulation strategy

This paper, as a part of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2b), assesses the impacts of different levels of global warming on the thermal structure of Lake Erken (Sweden). The General Ocean Turbulence Model (GOTM) one-dimensional hydrodynamic model was used to simulate water temperature when using ISIMIP2b bias-corrected climate model projections as input. These projections have a daily time step, while lake model simulations are often forced at hourly or shorter time steps. Therefore, it was necessary to first test the ability of GOTM to simulate Lake Erken water temperature using daily vs hourly meteorological forcing data. In order to do this, three data sets were used to force the model as follows: (1) hourly measured data, (2) daily average data derived from the first data set, and (3) synthetic hourly data created from the daily data set using generalised regression artificial neural network methods. This last data set is developed using a method that could also be applied to the daily time step ISIMIP scenarios to obtain hourly model input if needed. The lake model was shown to accurately simulate Lake Erken water temperature when forced with either daily or synthetic hourly data. Long-term simulations forced with daily or synthetic hourly meteorological data suggest that by the late 21st century the lake will undergo clear changes in thermal structure. For the representative concentration pathway (RCP) scenario, namely RCP2.6, surface water temperature was projected to increase by 1.79 and 1.36 ^∘C when the lake model was forced at daily and hourly resolutions respectively, and for RCP6.0 these increases were projected to be 3.08 and 2.31 ^∘C. Changes in lake stability were projected to increase, and the stratification duration was projected to be longer by 13 and 11 d under RCP2.6 scenario and 22 and 18 d under RCP6.0 scenario for daily and hourly resolutions. Model changes in thermal indices were very similar when using either the daily or synthetic hourly forcing, suggesting that the original ISIMIP climate model projections at a daily time step can be sufficient for the purpose of simulating lake water temperature.

1 Introduction

The thermal structure of lakes is controlled by heat and energy exchange across the air–water interface, which is in turn determined by meteorological forcing (Woolway et al., 2017). Climate change will affect air–water energy exchanges and alter the temperature regime and mixing of lakes (Woolway and Merchant, 2019). For example, increases in air temperature result in a consequent warming of lake water temperature (Sahoo et al., 2015) causing shorter ice-cover periods (Kainz et al., 2017; Butcher et al., 2015), longer stratified periods (Ficker et al., 2017; Woolway et al., 2017; Magee and Wu, 2017), and increased lake stability (Rempfer et al., 2010; Hadley et al., 2014). Decreasing wind speed can induce more stable and long-lasting stratification and increased epilimnetic temperature (Woolway et al., 2017, 2019).

The most direct effect of climate change on lakes is a warming of the lake surface temperature. For example, global average warming rates of 0.34 ^∘C per decade have been observed between 1985 and 2009 by O'Reilly et al. (2015). Hypolimnetic temperature responds less clearly to warming and has been observed to be warming, cooling or not changing significantly with increasing air temperature (Shimoda et al., 2011; Butcher et al., 2015; Winslow et al., 2017). And, these changing water temperatures have also led to an increased stability and duration of stratification (Butcher et al., 2015; Kraemer et al., 2015). A final consequence of warming lake temperature is projected to be the shift in the mixing regime (Kirillin, 2010; Shimoda et al., 2011; Shatwell et al., 2019; Woolway and Merchant, 2019). For example, loss of ice cover in deep lakes is likely to turn amictic lakes into cold monomictic lakes and cold monomictic lakes into dimictic lakes (Nõges et al., 2009). These changes in lake water temperature and thermal stratification influence lake ecosystem dynamics (MacKay et al., 2009).

Numerical modelling plays a key role in estimating the sensitivity of the lakes to changes in the climate. One-dimensional lake models are widely used due to their computational efficiency and the realistic temperature profiles they produce. Several studies have investigated the impacts of climate change on lake water temperature under the regional climatic model or global climatic model (RCM or GCM) projections (Persson et al., 2005; Kirillin, 2010; Perroud and Goyette, 2010; Samal et al., 2012; Ladwig et al., 2018; Shatwell et al., 2019; Woolway and Merchant, 2019). Commonly, when undertaking climate change impact studies, hydrodynamic lake models are driven by daily resolution RCM or GCM outputs. Bruce et al. (2018) undertook a comparative analysis of model performance, using daily and hourly resolution meteorological forcing data, and found a better agreement between observations and predictions of full-profile temperature when the lakes were modelled using hourly meteorological input. This reinforces the importance of diurnal forcing on one-dimensional model predictive capability.

The purpose of this study is therefore (1) to test the ability of a one-dimensional hydrodynamic model (General Ocean Turbulence Model – GOTM) to simulate the water temperature of Lake Erken (Sweden) using daily vs hourly meteorological forcing data for the period 2006–2016; (2) to develop a reliable method to disaggregate daily meteorological data to an hourly synthetic product that can be used to force the GOTM model and convert the daily GCM outputs available from the ISIMIP into hourly meteorological data sets; and (3) to assess the impacts on the thermal structure of Lake Erken at different levels of global warming when GOTM is driven by hourly and daily model projections. In fulfilling these objectives, this study provides the first evaluation of modelling methods that will be used by the lake sector within the ISIMIP.

2 Material and Methods

2.1 Study site

Lake Erken (59^∘51^′ N, 18^∘36^′ E) is a mesotrophic lake located in eastern central Sweden, with a maximum depth of 21 m, a mean depth of 9 m and a surface area of 23.7 km². The lake is dimictic with summer stratification usually beginning in May–June and ending in August–September, while the onset of ice cover occurs between December and February and ice loss is in April–May (Persson and Jones, 2008). It is the lake's relatively shallow depth and large surface area which leads to large interannual variability in the timing and patterns of thermal stratification since heat can be readily transferred through the shallow water column by wind mixing (Magee and Wu, 2017), and since the lake has a relatively low heat storage and, therefore, responds more directly to short-term variations in weather. The lake has a retention time of approximately 7 years and shows annual variations in water level that are less the 1 m (Pierson et al., 1992; Moras et al., 2019).

2.2 Lake model

General Ocean Turbulence Model (GOTM) is a one-dimensional water column model that simulates the most important hydrodynamic and thermodynamics processes related to vertical mixing in natural waters (Umlauf and Burchard, 2005). GOTM was developed by Burchard et al. (1999) for modelling turbulence in the oceans, but it has been recently adapted for use in hydrodynamic modelling of lakes (Sachse et al., 2014). The strength of GOTM is the vast number of well-tested turbulence models that have been implemented, spanning from simple prescribed expressions for the turbulent diffusivities to complex Reynolds stress models with several differential transport equations. Typically GOTM is used as a stand-alone model for investigating the dynamics of boundary layers in natural waters, but it can also be coupled to a biogeochemical model using the Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014).

2.3 Data sets

Local meteorological variables were collected either from a small island 500 m offshore from the Erken Laboratory or the Swedish Meteorological Hydrological Institute (SMHI) Svanberga station just behind the laboratory. The Malma Island meteorological station (59.83909^∘ N, 18.629558^∘ E) measured air temperature at 2 m above water surface, wind speed at 10 m above the water surface and short-wave radiation. These data were measured at 1 min intervals and saved as 60 min mean values. Mean sea level, pressure, relative humidity and precipitation were measured at the Svanberga meteorological station at 800 m from the Malma Island meteorological station (59.8321^∘ N, 18.6348^∘ E) with a frequency of 60 min. Hourly cloud cover was recorded from Svenska Högarna station (59.4445^∘ N, 19.5059^∘ E) at 69 km south-east of Lake Erken.

The measured hourly meteorological data were used to construct two additional data sets that would replicate the data resolution that could potentially be used to force the GOTM model with ISIMIP scenarios. First, to test running the model at a daily resolution, a daily data set was created by averaging the hourly one (except for precipitation which was summed). Second, this mean daily data set was disaggregated to form a synthetic hourly data set. Hourly estimations of air temperature, wind speed, relative humidity and short-wave radiation were estimated using the generalised regression artificial neural network (GRNN) methods described below. For atmospheric pressure and cloud cover, mean daily values were assumed to be constant over the day. Precipitation was disaggregated, assuming a uniform distribution of the daily total (Waichler and Wigmosta, 2003).

Since both of these data sets are based on the same measured hourly data, a comparison of model simulations of lake water temperature allows the importance of hourly vs daily temporal resolution in the forcing data to be evaluated, and also the improvements in model performance that can be obtained from daily data (as in the ISIMIP scenarios) when imposing a diurnal cycle on the mean daily data.

Water temperature data needed to calibrate the model was monitored from an automated floating station (59.84297^∘ N, 18.635433^∘ E). During ice-free conditions, measurements were made every 0.5 m from 0.5 m to a depth of 15 m. Measurements were made every minute, and a mean of these measurements was stored every 30 min.

2.4 Climate scenarios

The ISIMIP climate scenarios are bias-corrected global climate model (GCM; Hempel et al., 2013) data made available at daily temporal and 0.5^∘ horizontal resolutions for the variables listed in Table 1. All data needed as input to the GOTM model are available in these climate scenarios with the exception of cloud cover, which was estimated from short-wave radiation (Martin and McCutcheon, 1999). Data from the grid box overlying Lake Erken were available from the GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR and MIROC5 GCM models that were each run under three emission scenarios. These included a scenario with historical levels of atmospheric CO₂ between 1861 and 2005 and two future scenarios (RCP2.6 and RCP6.0) from 2006 to 2100. RCP2.6 is the strongest mitigation pathway that is expected to limit mean global warming to between 1.5 and 2 ^∘C. RCP6.0 is an intermediate mitigation pathway where global warming is projected to rise to between 2.5 and 4 ^∘C by the end of the century compared to the pre-industrial period (Frieler et al., 2017).

Table 1 Bias-corrected variables in the ISIMIP data set.

2.5 Temporal disaggregation of daily meteorological forcing data

Kathib and Elmernreich (2015) proposed a generalised regression artificial neural network (GRNN) model for predicting hourly variations in short-wave radiation from daily average measurements. Using the GRNN model to predict hourly solar radiation required 10 geographical and climatic variables as input including hour, day, month, latitude, longitude, daily average short-wave radiation, daily precipitation, the solar elevation associated with the hour, and time of sunrise and sunset. Precipitation was used to define wet and dry status that affected atmospheric attenuation (Waichler and Wigmosta, 2003).

There are also empirical models developed for calculating hourly air temperature, wind speed and relative humidity. Parton and Logan (1981) proposed a model for predicting diurnal variations in air temperature. Daylight air temperature was modelled using a sine wave with the minimum value at sunrise, maximum value at solar noon and mean value at sunset. Night-time air temperature was modelled as a linear interpolation between air temperature of the previous day and sunrise air temperature of the following day. Guo et al. (2016) generated hourly values of wind speed by computing a cosine function dependent on the mean daily wind speed, the maximum daily wind speed and the hour of the day when the wind speed is maximum. Waichler and Wigmosta (2003) estimated hourly values of relative humidity from daily maximum and minimum air temperature and daily maximum and minimum relative humidity. Using these studies as a guide, we developed GRNN models to predict hourly air temperature, wind speed and relative humidity. The input parameters for GRNN models were geographical variables and meteorological variables. The geographical variables include longitude, latitude, solar elevation associated with the hour, time of sunrise and sunset, hour, day and month, while the meteorological variables include average, maximum and minimum daily air temperature, daily wind speed, daily relative humidity, and daily precipitation.

The GRNN models were constructed using 8 years of data. From this whole set of data, the first 5 years, from 2008 to 2012, were used for training, and the final 3 years of data, from 2013 to 2015, were used for validating the results. The accuracy of the trained network was assessed by comparing the simulated output with actual observed hourly data. The performance index for training and validating sets of GRNN models is given in terms of mean bias error (MBE), root mean squared error (RMSE) and Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970). More detailed descriptions of the GRNN methods and models are given in Sect. S1 of the Supplement. The GRNN models were used to disaggregate the mean daily measured data, used to evaluate the necessity of disaggregation (Sect. 3.2), and also for all GCM scenarios (Sect. 3.3) to further evaluate the effects of disaggregation on the results of simulations of future changes in lake thermal structure.

2.6 Model set-up, calibration and validation

The GOTM model version 5.1 was used in this study. The meteorological parameters for running the model were air temperature (^∘C), wind speed (m s⁻¹), short-wave radiation (W m⁻²), cloud cover (dimensionless, 0–1), relative humidity (%), atmospheric pressure (hPa), and precipitation (mm d⁻¹ or mm h⁻¹). Inflows and outflows were not included in this study, and water level was considered fixed in the simulations. This version of GOTM did not have the ability to simulate lake ice, so for this study the inverse stratification period was not analysed. Moras et al. (2019) have shown that, despite this limitation, the model is able to accurately simulate water temperature and the phenology of thermal stratification during the remainder of the year. The initial conditions for water temperature were derived from a measured vertical profile. GOTM was run at an hourly model computational time step, and simulated water temperature was saved as daily mean values each 0.5 m (42 layers).

Calibration of the GOTM model was conducted to adjust the model parameters within their feasible range in order to minimise the error between measured and modelled temperature (Huang and Liu, 2010). A period of 9 years was selected for the calibration of GOTM, namely 2006–2014 (including 1 year spin-up time followed by 8 years for calibration). The model parameters that were calibrated were surface heat flux factor (shf_factor), short-wave radiation factor (swr_factor), wind factor (wind_factor), minimum turbulent kinetic energy (k_min), and e-folding depth for visible fraction of light (g2). The programme used to calibrate the model was Auto-Calibration Python (ACPy), developed by Bolding and Bruggeman (https://bolding-bruggeman.com/portfolio/acpy/, last access: 12 June 2018); it uses a differential evolution algorithm which calculates a log likelihood function based on comparing the modelled and measured water temperature (Storn and Price, 1997). The validation period was 2 years (2015–2016).

Table 2 GOTM lake model parameters and calibrated values.

For both calibration and validation, daily average water temperatures were simulated when GOTM was forced using the three meteorological data sets described above, namely measured average daily, measured average hourly and synthetic hourly data. Model simulated profiles of mean daily water temperature were then compared to mean daily measured water temperature. Three separate model calibrations were made based on simulations forced with the different meteorological data sets. During calibration the model was run approximately 10 000 times to obtain a stable solution specifying the optimal parameter set. The details of the feasible range of model parameters and the parameters associated with each calibration are given in Table 2.

Model performance was evaluated by comparing average daily modelled and measured temperature profiles and other metrics describing the lake thermal structure (surface and bottom temperature, volumetrically weighted averaged whole-lake temperature, Schmidt stability, thermocline depth, duration, and onset and loss of thermal stratification). The coefficients used to quantify the strength of model fit were MBE, RMSE and NSE.

2.7 Thermal indices

A range of thermal metrics, namely surface temperature (shallowest observation), bottom temperature (deepest observation) and thermocline depth (depth of the maximum density gradient), were derived on a daily basis from the daily simulated lake temperature profiles (temperature data with a vertical resolution of 0.5 m from 0.5 to 15 m depth). Also, Schmidt stability (resistance to mechanical mixing due to the potential energy inherent in the density stratification of the water column; Schmidt, 1928; Idso, 1973) and whole-lake temperature (volumetric weighted mean whole-lake temperature) were estimated from these profiles using Lake Analyzer (Read et al., 2011). The duration of thermal stratification was calculated as the longest continuous period when the water column density difference from the bottom to the surface of the lake was greater than 0.1 kg m⁻³ (according to ISIMIP2b lake sector protocol). The date of the onset and loss of the thermal stratification was defined as the first time that this density difference persisted for more than 5 d or was absent for at least 5 d (Kraemer et al., 2015).

2.8 Statistical analysis

Anomalies were calculated to further evaluate the impacts on lake water temperature and thermal stratification. The anomalies were computed for each GCM by taking the difference between the annual average of each year (2011–2100) from RCP2.6 and 6.0 scenarios and the average for the entire period 1981–2010 from the historical scenario. These average values were calculated using the months between April and September. The slope of the significant trends was evaluated by least squares linear regression, except when the residuals did not follow a normal distribution. Then the non-parametric Mann–Kendall test for the significance of trends and the Theil–Sen method (Theil, 1950; Sen, 1968) to estimate the slope of the significant trends were used instead. The Student t mean difference test was used to compare average values of each of the thermal indices. Distribution normality and variance homoscedasticity were assessed by the Shapiro–Wilk test and Fisher's F test respectively. When thermal indices time series did not follow a normal distribution the non-parametric Mann–Whitney U test (equal variances) or Kolmogorov–Smirnov test (different variances) were used instead. The statistical analysis was carried out using R version 3.4.4. The progress of climate-related impacts on the thermal stratification of the lake over the 21st century was assessed as the difference in mean conditions between the reference period (1981–2010) and both the mid-century (2041–2070) and the late century (2071–2100). Climate model data followed the same statistical analysis.

3 Results

3.1 Hourly meteorological modelling

There was a close agreement between GRNN model predictions and measured meteorological data as shown in Fig. 1, for 1 year, and in Sect. S1. For air temperature, short-wave radiation and humidity, the statistics of model performance always suggested a strong model fit in the training data sets and also remained strong, but somewhat lower, in the validation data sets (Table 3). NSEs were 1.00 for the training data sets and ranged from 0.69 to 0.94 for the validation data sets. Estimates of bias were very small. Wind speed was the variable showing the poorest performance with a NSE of 0.78 and 0.58 and RMSE values of 1.06 and 1.37 m s⁻¹ for the training and validate data sets. In general, the GRNN model tended to overestimate wind speed (MBE of 0.63±0.92 m s⁻¹) when the observed wind speed was lower than or equal to 3.84 m s⁻¹, whereas projected wind speed tended to be underestimated (MBE of $-\mathrm{0.78}\pm \mathrm{1.17}$ m s⁻¹) when the observations were greater than 3.84 m s⁻¹.

Figure 1 GRNN temporal disaggregation of meteorological forcing data. Observations vs simulations of (a) air temperature, (b) short-wave radiation, (c) relative humidity, and (d) wind speed for 2015 (validation data set).

Table 3 GRNN models' performance evaluation.

3.2 Lake model performance

Temperature observations and simulations for the three configurations of meteorological forcing data for both calibration and validation periods are shown in Fig. 2 and Sect. S2. Model performance metrics associated with these simulations are provided in Table 4. These data demonstrate that GOTM was able to accurately reproduce the measured temperature profiles. For an average of all three calibration data sets, a RMSE of 0.95 ^∘C and NSE of 0.94 were obtained. Temperature simulations in the shorter and less variable validation period (RMSE of 0.66 ^∘C and NSE of 0.97) were more accurate than for the calibration period, but in both periods the model performance was considered strong. For a full-profile temperature, the maximum RMSE value was 1.04 ^∘C and the minimum NSE was 0.93. Bottom temperature was least accurately simulated, with RMSE and NSE values reaching 1.33 ^∘C and 0.83 respectively.

Figure 2 GOTM water temperature simulations. Daily averaged water temperature in Lake Erken for the following calibration (1a–1d) and validation (2a–2d) periods: observations (1a, 2a), simulations driven by daily meteorological data (1b, 2b), hourly meteorological data (1c, 2c), and synthetic hourly meteorological data (1d, 2d).

Table 4 GOTM lake model performance evaluation of MBE, RMSE, and NSE for full-profile temperature, surface temperature, bottom temperature, volumetrically weighted averaged whole-lake temperatures, Schmidt stability, thermocline depth, duration, and onset and loss of thermal stratification using simulated results from running GOTM driven by daily (24 h met), hourly (1 h met), and synthetic hourly (synthetic 1 h met) meteorological data sets.

When comparing the metrics of model fit associated with simulations forced with the three different input data sets the simulations forced with mean daily input were slightly less accurate than those forced with either the measured or synthetic hourly input. As an example, full-profile RMSE values for the calibration period ranged from 0.88 to 1.04 ^∘C, with the lower error levels associated with simulations driven by hourly meteorological data sets, whereas for the validation period the RMSEs were comparable for all data sets. The MBE values of the full-profile temperatures indicated a slight cold temperature bias (average MBE of −0.05 ^∘C). The model performance predicting just surface temperatures was similar for all of the three calibrations (average RMSE of 0.67 ^∘C and NSE of 0.97) and was more accurate than the estimations of the full-profile temperatures. The MBE showed that for surface temperature GOTM also tended to produce a small cold temperature bias (average MBE of −0.10 ^∘C). The simulations of bottom temperature were slightly less accurate, having average RMSE of 0.96 ^∘C and NSE of 0.90, and also showed a tendency have a slightly higher RMSE values for calibrations forced with daily input. Also, the bottom temperature showed lower RMSE values for the validation period (average RMSE of 0.67 ^∘C) than the calibration period (average RMSE of 1.25 ^∘C), but, in contrast to the surface temperature, there was a slight warm temperature bias (average MBE of 0.06 ^∘C). When evaluating the simulations of volumetrically weighted averaged whole-lake temperatures, we found that model errors were of a similar magnitude for all simulations in both the calibration and validation periods with an average RMSE of 0.53 ^∘C and NSE of 0.98, tending to a slight cold temperature bias (average MBE of −0.08 ^∘C).

The calculation of Schmidt stability was also well simulated using all three data sets (average RMSE of 17.24 J m⁻² and NSE of 0.88). The lowest RMSE values were for the validation period (average RMSE of 13.34 J m⁻²), whereas during the calibration period the values were slightly greater (average RMSE of 21.14 J m⁻²). Thermocline depth was the parameter with the poorest performance (average RMSE of 3 m). The MBE values (average MBE of 0.80 m) indicate a bias towards an underprediction of thermocline depth (shallower thermocline depths). The RMSE associated with the prediction of the duration of stratification was, on average, 10.43 d, with lower RMSE values for the validation period (average RMSE of 8.04 d) than the calibration period (average RMSE 12.81 d). The simulations of the onset of the stratification were more accurate, having average RMSE of 2.64 d, but, in contrast, predictions of the loss of stratification were less accurate (average RMSE of 7.99 d).

3.3 Climate data projections

The lake model simulations undertaken here were forced by four climate model projections (GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR and MIROC5) that were in turn forced with three emissions scenarios (historical, RCP2.6 and RCP6.0). Average annual air temperature of the climate model ensemble for the reference period (1981–2010) was 11.88 ^∘C. Disaggregation of the climate input to an hourly time step resulted in a slightly warmer temperature (12.05 ^∘C)¹. Under future-scenario RCP2.6, the average increase was projected to be 2.22 ^∘C (1.71 ^∘C) by the mid-century (2041–2070), with a negligible change after the mid-century, as would be expected from this scenario with the strongest mitigation. During the period 2011–2100, air temperature increased at a rate from 0.08 to 0.17 ^∘C per decade (0.06 to 0.14 ^∘C per decade). In contrast, under RCP6.0, the average air temperature increased by 2.61 ^∘C (2.01 ^∘C) by the mid-century and continued rising to 3.61 ^∘C (2.76 ^∘C) by the late century. For RCP6.0 the trend in air temperature increased over the entire 2011–2100 period, on average, by 0.34 ^∘C per decade (0.26 ^∘C per decade) over all GCMs, with the individual trends ranging from 0.18 to 0.43 ^∘C per decade (0.14 to 0.33 ^∘C per decade). For the remaining meteorological variables there were fewer distinct changes between the historical and future periods. Under the RCP2.6 scenario the overall annual mean change in wind speed was negligible, while under RCP6.0 two options were projected, namely an increase (GFDL–ESM2M and MIROC5) and a decrease (HadGEM2–ES and IPSL–CM5A–LR). Relative humidity was projected to decrease for future-scenario RCP2.6 and 6.0. For RCP6.0, significant trends ranged from 0.29 % to 0.36 % per decade in the interval 2011–2100. An increase in short-wave radiation was projected for all RCP scenarios by the late century, with a negligible mean change after the mid-century under RCP6.0. The increase in short-wave radiation is coupled with a decrease in cloud cover. By the late century, the mean decrease in cloud cover was projected to be 0.06 for RCP2.6 and 0.07 for RCP6.0. More detailed evaluations of the differences in the climate projection based on the original ISIMIP daily time step and the hourly disaggregated data are given in Sect. S3.

3.4 Long-term modelled changes in thermal stratification

Lake model simulations were made using both the original daily resolution of the ISIMIP–GCM scenarios and also at hourly resolution using disaggregated data developed using the GRNN models. Simulated water temperatures for the historical, RCP2.6 and 6.0 scenarios under daily IPSL–CM5A–LR projections are presented as temperature isopleths in Fig. 3 and Sect. S4. These were created by averaging the daily temperature profiles for all years associated with each of the emission scenarios. These mean scenario temperature isopleths provide a clear visualisation of how, for future scenarios, surface and bottom water temperatures are projected to increase with stronger and shallower stratification, an earlier stratification onset, a later fall overturn and, consequently, a longer stratification period.

Figure 3 Temperature isopleth diagrams for the (a) historical, (b) RCP2.6, and (c) RCP6.0 scenarios and showing results from the lake model forced with daily IPSL–CM5A–LR projections. The temperature matrix used to make these plots was created by averaging the simulated daily temperature profiles for every year in each scenario.

Figure 4 Evolution of annual average projected anomalies (from April to September) for (1a, 2a) whole-lake temperature, (1b, 2b) surface temperature, (1c, 2c) bottom temperature, (1d, 2d) Schmidt stability, and (1e, 2e) thermocline depth and showing results from when the lake model was forced with daily GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR, and MIROC5 projections from 2011 to 2100 under RCP2.6 and 6.0. Anomalies are relative to the reference period (1981–2010).

Figure 5 Evolution of annual average projected anomalies (from April to September) for (1a, 2a) duration and (1b, 2b) onset and (1c, 2c) loss of stratification and showing results from when the lake model was forced with daily GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR, and MIROC5 projections from 2011 to 2100 under RCP2.6 and 6.0. Anomalies are relative to the reference period (1981–2010).

Table 5 Trend analysis from 2011–2100 for surface temperature, bottom temperature, whole-lake temperature, Schmidt stability, thermocline depth, duration, and onset an d loss of stratification (ns: not significant) for RCP6.0. Italics denote trends with p-value > 0.05.

In Figs. 4–5 we show the long-term trends in the anomalies in lake thermal metrics simulated to occur over the RCP2.6 and 6.0 emission scenarios. Trends in whole-lake temperature calculated for over a period of 90 years (2011–2100) were projected to increase except in the case of GFDL–ESM2M, which showed weaker or non-significant changes for all measures of thermal stratification (Table 5 and Sect. S4). Under RCP2.6, significant trends ranged from 0.07 to 0.10 ^∘C per decade (0.05 to 0.08 ^∘C per decade), but most of the change occurred in the first half of the century. For RCP6.0, the projected rate of change ranged from 0.14 to 0.26 ^∘C per decade (0.10 to 0.19 ^∘C per decade). By the late century, the mean projected increase in whole-lake temperature was 1.34 ^∘C (1.00 ^∘C) for RCP2.6 and 2.39 ^∘C (1.75 ^∘C) for RCP6.0, with a negligible change after the mid-century under RCP2.6 (Fig. 6, Table 6 and Sect. S4).

Figure 6 Changes in annually averaged thermal metrics (from April to September) for (2a, 3a) whole-lake temperature, (2b, 3b) surface temperature, (2c, 3c) bottom temperature, (2d, 3d) Schmidt stability, and (2e, 3e) thermocline depth under RCP6.0 and showing results from when the lake model was forced with daily GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR, and MIROC5 projections. All changes for the mid-century (2041–2070) and the late century (2071–2100) are relative to the reference period (1981–2011). The mean (vertical line) is also shown. Changes in thermal metrics greater than 0 show an increase and lower than 0 show a decrease.

Table 6 Average thermal metrics for the reference period (1981–2010), and average projected change in thermal metrics for the mid-century and the late century for RCP6.0.

Changes in lake surface temperature were, as expected, greater than for whole-lake temperature. For the reference period, the mean April–September surface temperature was on average 13.61 ^∘C (13.84 ^∘C) warming up significantly over the 21st century, so that by the late century the average projected increase was 1.79 ^∘C (1.35 ^∘C) for RCP2.6 and 3.08 ^∘C (2.31 ^∘C) for RCP6.0. From 2011 to 2100 there was a significant long-term trend for RCP2.6 surface temperature, which increased at a rate from 0.07 to 0.15 ^∘C per decade (0.06 to 0.13 ^∘C per decade). Under RCP6.0 the mean surface temperature increase of the ensemble was 0.30 ^∘C per decade (0.23 ^∘C per decade) ranging between 0.16 and 0.38 ^∘C per decade (0.12 to 0.29 ^∘C per decade). The projected increase in bottom temperature was not as marked as it was for the other metrics of lake temperature. On average, the bottom temperature increased from 9.20 ^∘C (9.67 ^∘C) in the reference period to 9.77 ^∘C (9.99 ^∘C) and 10.32 ^∘C (10.34 ^∘C) by the late century for RCP2.6 and 6.0 respectively. Significant rates of change in bottom temperature were not predicted during the RCP2.6 scenario, but for the RCP6.0 scenario bottom temperature did undergo significant warming rates for HadGEM2–ES and MIROC5 projections were 0.06 ^∘C per decade and 0.11 ^∘C per decade (0.09 ^∘C per decade) respectively.

Figure 7 Changes in annually averaged thermal metrics (from April to September) for (2a, 3a) duration and (2b, 3b) onset and (2c, 3c) loss of stratification under RCP6.0 and showing results from when the lake model was forced with daily GFDL–ESM2M, HadGEM2–ES, IPSL–CM5A–LR, and MIROC5 projections. All changes for the mid-century (2041–2070) and the late century (2071–2100) are relative to the reference period (1981–2011). The mean (vertical line) is also shown. Changes in thermal metrics greater than 0 show an increase and lower than 0 show a decrease.

There were also projected changes in the resistance to mechanical mixing. For the reference period, an average of 68.65 J m⁻² (65.56 J m⁻²) was required to completely mix the water column, while by the late century it increased by 29.08 J m⁻² (22.74 J m⁻²) for RCP2.6 and 49.22 J m⁻² (38.07 J m⁻²) for RCP6.0 (Fig. 4, Table 6 and Sect. S4). This greater stability also corresponds to a longer duration of stratification. From 2011 to 2100, a significant increase in the duration stratification was projected for both future scenarios RCP2.6 and 6.0, ranging from 1.13 to 1.70 d per decade (0.87 to 1.30 d per decade) for RCP2.6 and 2.45 to 3.56 d per decade (2.00 to 3.08 d per decade) for RCP6.0 (Fig. 5, Table 5 and Sect. S4), which led to a mean change of 13 d (11 d) and 22 d (18 d) for RCP2.6 and 6.0 respectively (Fig. 7, Table 6 and Sect. S4). The longer period of stratification resulted from both an earlier onset of thermal stratification and a later loss of thermal stratification (Figs. 5 and 7, Tables 5–6 and Sect. S4). Mean annual thermocline depth was simulated to be shallower under future conditions. By the late century, the reduction in thermocline depth was projected to be 0.38 m (0.41 m) for RCP2.6 and 0.49 m (0.57 m) for RCP6.0, although a significant trend in the decline was only found for the later scenario.

The trends in Figs. 4–5 are quite variable from year to year, and as would be expected, there is no direct correspondence in the temporal variations of one GCM relative to another. To provide an alternative method of comparing the changes simulated by the future climate scenarios shown in Figs. 4–5, the daily anomalies for each trend line are also presented as frequency distributions in the Figs. 6–7 for the simulations made under the RPC6.0 scenario. These show that for all metrics there is a clear shift in the lake thermal conditions that are consistent with a warmer climate, but also that there are extremes in the distributions that can lead to unrepresentative results when, for example, future conditions briefly return to historical levels, or when the effects of warming are much greater that would be expected on average. This later case can cause important changes in lake ecology if the extreme conditions result in a change in lake state by the passing of a tipping point. Figures 6–7 also clearly show the differences in simulations forced by different GCMs. Most obvious is the difference in the results from GFDL–ESM2M, which consistently simulated smaller changes in lake thermal structure during the mid- and late-century periods despite having a data distribution that was similar to the other models during the historical period.

3.5 Comparison between long-term thermal metrics derived from daily and hourly climate data

Future changes in thermal metrics based on both RCP2.6 and RCP6.0 were slightly greater when the GOTM model was forced at daily resolutions (Tables 5–6 and Sect. S4) than at an hourly resolution. This included changes in mean surface temperatures and also in the annual average whole-lake temperature (Sect. S5). However, under RCP2.6 non-significant differences were found for bottom temperature, Schmidt stability, thermocline depth, the duration, and onset and loss of stratification. In all cases where differences were found between the simulations forced with daily vs hourly data there were no changes in direction and only minor changes in the magnitude of the change suggested by the simulations (Sects. S4–S5).

4 Discussion

The simulated water temperature and related metrics of thermal stratification were in excellent agreement with the extensive set of measured water temperature data that were available for model calibration at Lake Erken (Moras et al., 2019; Fig. 3; Table 3). Water temperature simulations were apparently more accurate for the validation period (2015–2016) than for the calibration period (2006–2014), which may appear unusual, but this is due to the higher variability in observed water temperature during the longer calibration period. Years with a longer duration of stratification and stronger stability, generally had higher simulation errors. Half of the 8-year calibration period exhibited these conditions, while the 2 years used for validation both exhibited shorter duration of stratification and weaker stability. The thermocline depth was the thermal metric that was predicted with the greatest uncertainty. This is in part caused by the presence of internal seiches in Lake Erken, which result in the measured temperatures in the region of the thermocline having a level of variability that cannot be reproduced by one-dimensional models such as GOTM. Bruce et al. (2018) detected a strong correlation between the accuracy of the extinction coefficient and model simulations of full-profile temperature and thermocline depth, which shows the importance of light extinction in the prediction of thermocline depth. Since a single fixed value of e-folding depth (Table 2) for the visible fraction of the light (the inverse of the extinction coefficient) was used in the GOTM simulations, the effects of seasonal variations in light extinction (Perroud et al., 2009) on thermocline depth were not evaluated.

The model parameters adjusted during the calibration processes were non-dimensional scaling factors (shf_factor, swr_factor and wind_factor) and physical parameters which strongly influence the vertical distribution of light and temperature (k_min and g2). Wind is the dominant driver of mixing in lakes, and increases or decreases in wind speed (wind_factor) change the amount of turbulent kinetic energy available for mixing. The wind-scaling factor is often important when wind measurements occur some distance from the lake and/or accounts for wind-sheltering effects (Markfort et al., 2010). One would not expect the wind factor to deviate greatly from 1.0 at Lake Erken where wind is measured on an island in the lake. However, the dominant wind direction is along the lake's longest east–west fetch (Yang et al., 2016), which could explain the need to scale up the unidirectional wind speed measurements that were used as an input to GOTM. Furthermore, since it is the cube of wind speed that affects lake mixing, using longer averaging periods will underestimate the effects of gusting and variable winds. This could explain why we obtain higher calibrated values of the wind_factor when forcing the model with measured daily, rather than hourly, data (Table 2). Higher values of the wind_factor were also obtained when GOTM was forcing with synthetic hourly meteorological drivers. This is due to an underestimation of the faster wind speed predictions from the GRNN model (Fig. 1 and Sect. S1.). During the ACPy calibration, each of these parameters were calibrated while simultaneously influencing each other; shf_factor, swr_factor, wind_factor and g2 have a strong influence on heat and energy exchange across the air–water interface. There is to some extent an unavoidable tendency for the error in one parameter to be cancelled out by opposite errors in another parameter. This also illustrates how, to some extent, the calibration process can compensate for some of the limitations related to the temporal resolution of the model forcing data.

The performance of the GOTM model obtained in this study is comparable with results reported by others. Moras et al. (2019), who ran GOTM using hourly measured meteorology for a 57-year period, obtained a RMSE for daily full-profile water temperature of 1.09 ^∘C. Using the DYnamic REservoir Simulation Model (DYRESM), Magee and Wu (2017) reported RMSE values of 0.30 and 0.53 ^∘C for Lake Mendota and 1.45 and 1.94 ^∘C for Fish Lake for temperature estimates of the epilimnion and hypolimnion respectively. Perroud et al. (2009) simulated water temperature profiles of Lake Geneva over a 10-year period and obtained RMSE values of <2 ^∘C for the DYRESM model and <3 ^∘C for the Simstrat model.

The projected changes in lake thermal metrics depend on the selected GCM model and the future scenario or representative concentration pathway (RCP) that was simulated. The range of greenhouse gas (GHG) emissions included in this study was a stringent mitigation scenario (RCP2.6) and an intermediate scenario (RCP6.0). This is consistent with the ISIMIP2b simulation strategy that is intended to evaluate RCP2.6 as a scenario that aims to keep global warming below 2 ^∘C above pre-industrial temperatures by 2100. In contrast, for RCP6.0 increased levels of GHG emissions suggest that the global mean temperature will continually increase by 2.5 and 4 ^∘C by the end of the century. The effects of the mitigation measures that were adopted in RCP2.6 on lake thermal structure become most apparent in the late century. For example, for MIROC5 (when the lake model was forced at daily resolutions) the projected surface temperature change for the mid-century was similar for the two RCPs (2.10 ^∘C for RCP2.6 and 1.98 ^∘C for RCP6.0), but for the late-century period the projected change in surface temperature diverges among RCPs. Under RCP2.6 the surface temperature change declines from 2.10 to 1.80 ^∘C, while under RCP6.0 the change in surface temperature was projected to further increase from 1.98 to 2.97 ^∘C. Similar changes were projected for all thermal metrics. Under RCP6.0 there was also an increase in bottom temperature but at rates that were slower than surface temperature. Changes in lake stability increased from 38.67 J m⁻² by the mid-century to 64.62 J m⁻² by the late century, which increased the duration of stratification (from 16 to 22 d). While there was a general agreement among the models in the direction and relative magnitude of change in many of the metrics of lake thermal structure, there were also some differences among GCMs (Figs. 4–7 and Sect. S4), especially in relation to the GFDL–ESM2M model which consistently estimated lower levels of change. For example, by the late century the largest changes in surface temperature were projected for HadGEM2–ES (4.04 ^∘C) and the lowest were for GFDL–ESM2M (1.67 ^∘C) under future-scenario RCP6.0 when the lake model was forced at daily resolutions. However, the increase in the projected bottom temperature for GFDL–ESM2M (1.24 ^∘C) was greater than for HadGEM2–ES (0.91 ^∘C). This could be in part due to the projected changes in wind speed. The wind speed was projected to increase by 0.18 m s⁻¹ for GFDL–ESM2M, transferring heat to the lake bottom, but for HadGEM2–ES the wind speed decreased by 0.15 m s⁻¹ (atmospheric stilling; Woolway et al., 2017, 2019), reducing the magnitude of vertical mixing. This resulted in a greater gradient between surface and bottom temperatures and higher increases in the Schmidt stability (77.57 J m⁻²). This increased thermal gradient for HadGEM2–ES promoted shallower thermocline depth (1.26 m), but for GFDL–ESM2M a lower change in lake stability was projected (12.26 J m⁻²) leading to a deeper thermocline depth (0.22 m). Higher surface water temperature and stronger Schmidt stability can explain why the increased duration of stratification was projected to be longer for HadGEM2–ES (34 d) than for GFDL–ESM2M (6 d). The small change in thermal stability also explains why no change in loss of stratification was projected for GFDL–ESM2M. This illustrates the complexity of climate model–lake model interactions and clearly shows the importance of ensemble model simulations, as adopted by ISIMIP, for evaluating the effects of climate change on lakes.

When calibrating the GOTM model we found that model errors between simulated and measured water temperature were similar when GOTM was forced with either measured hourly or synthetic hourly meteorological data, and that the results obtained from the calibrations forced with mean daily metrological input were also similar to those obtained from the calibrations based on hourly input. This suggests that the daily time step of the ISIMIP climate scenarios is sufficient for forcing the GOTM model, and that for most studies within the ISIMIP lake sector disaggregation to hourly time steps will not be necessary. For example, changes in surface water temperature was to the order of 0.29 ^∘C per decade, with simulations forced with daily inputs, and 0.22 ^∘C per decade with hourly input data for MIROC5 under RCP6.0. These differences are of the same magnitude as the differences simulated using different GCM models. Similar levels of model performance using daily or hourly forcing data were obtained in part because of separate calibrations when the GOTM model was forced with the different data sets. Changes in the calibrated parameters used to characterise the lake thermal structure (Table 2) can apparently compensate for the lack of diurnal variability in the daily forcing data. GRNNs proved to be an effective method to disaggregate daily GCM forcing to an hourly temporal resolution for different weather variables such as air temperature, short-wave radiation, etc. However, GRNNs require a training phase in which the diurnal patterns to be learnt are presented to the network from historical meteorological measurements, and therefore if there are future changes in diurnal patterns, these cannot be reproduced. In addition, there is a high computational cost of disaggregating and storing the long-term daily climate data into an hourly data set.

The projected changes in thermal metrics were strongly influenced by the GMCs used to drive the water temperature simulations. Due to the high interannual variability, long periods of simulation were needed to ensure that the uncertainty is fully represented (Figs. 4–7; Sects. S3–S4). Under RCP6.0, trends in surface temperature calculated for the period 2011–2100 were projected to increase 0.38 ^∘C per decade for both HadGEM2–ES and IPSL–CM5A–LR when the lake model was forced at daily resolutions. However, 5th, 50th and 95th percentiles for surface temperature anomalies differ between models and are 0.84, 2.93, and 4.86 ^∘C for HadGEM2–ES and 0.33, 2.56, and 4.37 ^∘C for IPSL–CM5A–LR. Placing the probability density function (pdf) for HadGEM2–ES to the right of the pdf for IPSL–CM5A–LR illustrated that more extreme increases in surface temperature were projected by HadGEM2–ES. Projected bottom temperatures differed between HadGEM2–ES and IPSL–CM5A–LR. HadGEM2–ES was left-skewed and the median was 0.58 ^∘C, while IPSL–CM5A–LR pdf was right-skewed and the median was 1.16 ^∘C. As a consequence, lake stability was stronger for HadGEM2–ES (5th and 95th percentiles were 5.22 and 110.55 J m⁻²) than for IPSL–CM5A–LR (5th and 95th percentiles were −18.77 and 90.64 J m⁻²), even though the Schmidt stability medians were similar for both GCMs. A similar result occurred when projecting a longer duration of stratification for HadGEM2–ES (5th, 50th, and 95 percentiles were −0.63, 26.37, and 49.42 d) than IPSL–CM5A–LR (5th, 50th, and 95th percentiles were −10.33, 17.67, and 40.92 d). GCMs are useful for assessing climate change impacts on lakes, but GCM configurations vary from one to another. Therefore, it is crucial to assess different GCMs in a probabilistic manor (Figs. 4–7) to encapsulate the uncertainty in the lake thermal metrics without compromising the variability.

The study carried out by Moras et al. (2019) found changes in the phenology of Lake Erken's thermal stratification from 1961 to 2017. A significant increase in summer epilimnetic and whole-lake temperature of 0.35 and 0.24 ^∘C per decade occurred over the entire study period, while in spring and autumn larger significant positive trends were detected over the subinterval (1989–2017). In the present work future changes under the RCP6.0 emissions scenario found trends that were of a similar magnitude. The summer trends for the period 2011–2100 projected a significant increase in epilimnetic and whole-lake temperature ranging from 0.19 to 0.45 and 0.15 to 0.26 ^∘C per decade, respectively, when the lake model was forced at daily resolutions, while changes in summer hypolimnetic temperature were non-significant. During the spring and autumn significant increases in epilimnetic whole-lake temperature were also projected under RCP6.0 when the lake model was forced at daily resolutions, but they were somewhat lower than the ones detected by Moras et al. (2019). The increase in spring epilimnetic and whole-lake temperature ranged from 0.15 to 0.38 and 0.15 to 0.30 ^∘C per decade, while those in Moras et al. (2019) showed a higher rate of warming (0.44 and 0.40 ^∘C per decade for epilimnetic and whole-lake temperature respectively), and the GCM simulations promoted shorter durations in stratification. The projected increases in spring and autumn hypolimnetic temperatures were similar in magnitude, and in summer non-significant trends were detected both in this study and in Moras et al. (2019). The simulations made here and by Moras et al. (2019) are for the same lake and use the same lake model. The fact that the simulations presented here, using the RCP6.0 emission scenarios, showed similar or slightly lower rates of change compared to the simulations made by Moras et al. (2019) when the model was forced with measured historical data are unexpected given that the RCP6.0 scenario would project an accelerated rate of climate change compared to the historical period. This suggests that, at least for Lake Erken, future changes in lake thermal structure based on the ISIMIP2b–GCM projections may to some extent underestimate the actual changes that will occur.

The projected changes in thermal stratification can influence many aspects of the lake ecosystem. Increases in thermal stability and duration of stratification can intensify hypolimnetic oxygen depletion (Foley et al., 2012; Schwefel et al., 2016) and hence induce enhanced internal phosphorous loading (North et al., 2014), increase the release of dissolved iron and manganese from sediments (Schultze et al., 2017), and also increase methane emissions (Grasset et al., 2018). Warming lake temperature affects biological rates of metabolism, growth and reproduction (Rall et al., 2012) and can promote cyanobacterial blooms (Paerl and Paul, 2012). When coupled to a reduction in oxygen-rich water, warming water temperature leads to lower fish populations (O'Reilly et al., 2003; Yankova et al., 2017). Increase in evaporation associated with warming can lead to a decline in lake water levels (Hanrahan et al., 2010) with implications for water security.

5 Conclusion

In this study, which is the first test simulating lake hydrothermal structure following ISIMIP2b protocols, ensemble simulations show that changes in Lake Erken's surface temperature are projected to increase on average by 1.79 ^∘C for RCP2.6 and by 3.08 ^∘C for RCP6.0, and the length of the stratification is also projected to be longer by 13 d for RCP2.6 and by 22 d for RCP6.0 by the end of the 21st century. Most changes in the RCP2.6 scenario occurred during the first half of the century, suggesting that the aggressive mitigation methods represented in this scenario would be effective at reducing future changes in lake thermal structure. We also extensively document coinciding changes in water column temperatures, water column stability and thermocline depth both in this paper and in the Supplement. When combined, these results suggest important changes in the factors affecting lake biogeochemistry both directly, through changes in temperature, and indirectly, by influencing the availability of light and nutrients. By providing an initial test for the simulations that will be carried out by the ISIMIP lake sector, this paper sets the stage for more extensive worldwide evaluations of the effects of climate change on lake thermal structure.

This study showed the ability of the GOTM model to accurately simulate Lake Erken water temperature when the model was forced using either daily or hourly temporal resolution inputs. Neural networks were shown to be an effective method to disaggregate different weather variables, such as air temperature and short-wave radiation, in order to generate synthetic hourly meteorological data from the daily data that are typically available from GCM models. Model performance was slightly improved when using the synthetic hourly data, and climate change effects were somewhat lower when using such data to drive future climate simulations. However, it is not clear if data disaggregation is needed given the computational costs of developing such data sets and running long-term simulations at an hourly time step. Future climate simulations showed similar trends in the anomaly distributions when the GOTM model was forced with mean daily or synthetic hourly meteorological data, and we also found evidence that the calibration procedure partly compensates for differences in the temporal resolution of the model input.