A framework for ensemble modelling of climate change impacts on lakes worldwide: the ISIMIP Lake Sector

Empirical evidence demonstrates that lakes and reservoirs are warming across the globe. Consequently, there is an increased need to project future changes in lake thermal structure and resulting changes in lake biogeochemistry in order to plan for the likely impacts. Previous studies of the impacts of climate change on lakes have often relied on a single model forced with limited scenario-driven projections of future climate for a relatively small number of lakes. As a result, our understanding of the effects of climate change on lakes is fragmentary, based on scattered studies using different data sources and modelling protocols, and mainly focused on individual lakes or lake regions. This has precluded identification of the main impacts of climate change on lakes at global and regional scales and has likely contributed to the lack of lake water quality considerations in policy-relevant documents, such as the Assessment Reports of the Intergovernmental Panel on Climate Change (IPCC). Here, we describe a simulation protocol developed by the Lake Sector of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) for simulating climate change impacts on lakes using an ensemble of lake models and climate change scenarios for ISIMIP phases 2 and 3. The protocol prescribes lake simulations driven by climate forcing from gridded observations and different Earth system models under various representative greenhouse gas concentration pathways (RCPs), all consistently bias-corrected on a 0.5^∘ × 0.5^∘ global grid. In ISIMIP phase 2, 11 lake models were forced with these data to project the thermal structure of 62 well-studied lakes where data were available for calibration under historical conditions, and using uncalibrated models for 17 500 lakes defined for all global grid cells containing lakes. In ISIMIP phase 3, this approach was expanded to consider more lakes, more models, and more processes. The ISIMIP Lake Sector is the largest international effort to project future water temperature, thermal structure, and ice phenology of lakes at local and global scales and paves the way for future simulations of the impacts of climate change on water quality and biogeochemistry in lakes.

1 Introduction

There are over 117 million lakes on Earth covering only 3 % of the land surface (Verpoorter et al., 2014), yet freshwater ecosystems in general host 10 % of Earth's known animal species (Reid et al., 2019). Many lakes provide ecosystem services to their local communities including drinking water, fisheries, and transportation, and the number of services provided by lakes has been shown to decrease with deteriorating lake health (Janssen et al., 2021). In addition, lakes are effective as local indicators for both environmental changes at the watershed scale and as “sentinels of climate change” in that they buffer synoptic-scale variability but incorporate information on seasonal cycling, inter-annual variability, and long-term changes in lower atmospheric conditions. Therefore, studying lake impacts across scales is an important field of research for disentangling the global impacts of climate change from the other anthropogenic pressures that climate change interacts with. However, estimates of historical and future lake responses to climate change have, until recently, largely been carried out as site-specific studies with different goals, data, and modelling protocols, which complicates the generalization of simulated impacts at regional and global scales (Settele et al., 2015; IPCC, 2018).

Historical records show that lakes are already responding to climatic change by warming (O'Reilly et al., 2015; Pilla et al., 2020; Gal et al., 2020; Jane et al., 2021; Woolway et al., 2019b), experiencing declining ice cover (Weyhenmeyer et al., 2011; Sharma et al., 2019), shifting thermal habitats (Kraemer et al., 2021), changing mixing regimes (Kraemer et al., 2015; Woolway and Merchant, 2019), and decreasing oxygen levels (Jane et al., 2021). However, long-term monitoring data remain limited to a relatively small number of well-studied lakes, while time series from automated high sampling frequency monitoring buoys are still generally short (Marcé et al., 2016). The existing empirical evidence needs to be combined with lake models to understand how lakes have responded to historical changes (Moras et al., 2019) and how they could behave under future climatic change. Numerous numerical models have been used to assess climate change impacts on lake ecosystems (Schwefel et al., 2016; Hansen et al., 2017; Wang et al., 2018; Zwart et al., 2019; Ayala et al., 2020; Piccolroaz et al., 2020); however, the climate change impacts synthesized in the recent IPCC reports remain limited to a small number of lakes, regions, or specific impact models or climate change scenarios (Kraemer et al., 2015; IPCC, 2018). Multi-model ensemble simulations are increasingly used to obtain more robust assessments of freshwater ecosystem responses to climate change, but so far, only a few lakes have been assessed following a multi-model approach (Perroud et al., 2009; Trolle et al., 2014; Stepanenko et al., 2010, 2013, 2014; Thiery et al., 2014a; Gal et al., 2020; Guseva et al., 2020). To date, no multi-model ensembles have been used over a broad range of lakes to make either hindcast or future climate simulations, which would allow evaluation of the variability in model output related to different model formulations or parametrizations.

The ISIMIP framework (http://www.isimip.org, last access: 23 May 2022) provides a set of climate and socioeconomic forcing data to make consistent historical hindcast and future climate impact projections and evaluate impacts in response to policy-relevant climate change scenarios. ISIMIP is organized in different sectors ranging from hydrology to human health, all of which make use of common and openly provided input data. As part of ISIMIP, we initiated the Lake Sector and developed a lake model simulation protocol to assess climate change impacts on lakes and to provide robust scientific evidence of historical and potential future lake ecosystem changes. To this end, we used two complementary strategies: (i) a local strategy to simulate well-studied lakes where sufficient data were available for lake-specific model parameterization and calibration and (ii) a global strategy that applied lake models to simulate generic lakes for each lake-containing grid cell of the ISIMIP global grid. The simulation set-up described by the protocol enables impacts of climate change on lake characteristics to be projected and attributed. The protocol incorporates uncertainties from the differences in the global climate models (GCMs) providing forcing data, the differences in lake impact model structure, and lake geographical and ecosystem characteristics. The standardized output produced by the hydrodynamic lake models includes vertical profiles of water temperature and metrics describing thermal and ice conditions at daily to annual timescales. This multi-model ensemble provides a systematic overview of plausible future responses of lake ecosystems to a warming climate at an unprecedented geographical coverage. The forcing data from the GCM ensemble further enable the quantification of lake responses to changes in meteorological variables other than air temperature, e.g. from wind velocity or cloud cover (Woolway et al., 2019a), including their potential interactions with increases in air temperature. Completed and ongoing thermal regime simulations will provide the foundation for the modelling of water quality, greenhouse gas emissions, algal blooms, and water level fluctuations to be addressed in future ISIMIP rounds.

Here, we describe the protocol for the global- and local-scale intercomparison of lake model simulations completed for the second phase of ISIMIP (ISIMIP2), as well as the extensions to this protocol that have been implemented for the new phase three simulations (ISIMIP3). The evolution of the modelling protocol from ISIMIP2 to ISIMIP3, and the rationale for these advancements, will be described in individual sections related to the experimental set-up of the Lake Sector, such as changes to lake model forcing datasets and background information on lake mapping. First, we provide an overview of the climate data and climate change scenarios available through ISIMIP that were used as forcing data for lake impact models. Next, we explain the rationale behind the ISIMIP Lake Sector, give an overview of the impact simulations and briefly describe the lake models used at local and global spatial domains. Finally, we highlight some examples of the first lake impact model simulations from the ISIMIP2 simulation phase and illustrate how the simulations allowed us to quantify sources of uncertainty in future projections.

The Lake Sector simulations are the first ensemble projections using a consistent modelling framework to evaluate the impact of climate change on lakes, thereby informing researchers, policymakers, and water managers and enabling comparisons of impacts with other sectors participating in ISIMIP. Given the importance of lake thermal structure in regulating lake processes, ISIMIP2 simulation results from the Lake Sector provide a fundamental contribution to lake specific policy recommendations in reports from organizations like the IPCC, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES), and the UN Environment Programme World Water Quality Alliance (UNEP-WWQA). Furthermore, the work described here paves the way for lake water quality simulations playing a fundamental role in international policy.

2 The Lake Sector in the ISIMIP framework

The ISIMIP2 simulation framework was divided into two simulation rounds: ISIMIP2a and ISIMIP2b. ISIMIP2a focused on historical simulations, which were forced with gridded climate reanalysis products based on observed meteorological data. The ISIMIP2b simulation round focussed on quantifying the impacts of GCM-derived historical and CMIP5 (Coupled Model Intercomparison Project Phase 5) projected climate change relative to a pre-industrial control (Frieler et al., 2017). For the Lake Sector, ISIMIP2a simulations centred around the calibration of lake models at the local scale that were used for the future climate ensemble simulations in ISIMIP2b. These simulations also evaluated the lake models' ability to simulate observed climate variability and extremes at both local and global scales. The subsequent simulation phase, ISIMIP3a-b, will build on the latest gridded observations and CMIP6 global climate model simulations to provide meteorological forcing and improve the representation of non-climatic input data, such as land use and a dynamic global lake mask, which will be used to produce a new generation of lake model simulations.

The inception of the Lake Sector and the gathering of the first collection of models, modelling teams, and lake data providers greatly benefited from contributions from two global collaborative projects. First, the Lake Model Intercomparison Project (LakeMIP, https://www.unige.ch/climate/lakemip, last access: 23 May 2022) started in 2008 (Stepanenko et al., 2010) with the objective of comparing the thermodynamic regime of lakes (including lake-atmosphere interactions) in a wide range of climatic conditions and mixing regimes as simulated by several one-dimensional lake models (Perroud et al., 2009; Stepanenko et al., 2013, 2014; Thiery et al., 2014a; Guseva et al., 2020). Second, the Global Lake Ecological Observatory Network (GLEON, https://gleon.org, last access: 23 May 2022, Weathers et al., 2013) started in 2005, with the aim of sharing and interpreting lake data to understand, predict, and communicate the role and responses of lakes in a changing global environment.

3 Experimental set-up

The simulations followed the network-wide simulation protocols for ISIMIP 2a–b (Frieler et al., 2017; Schewe et al., 2019) and ISIMIP3a–b (see https://www.isimip.org/protocol/, last access: 23 May 2022 for an overview). Here, we describe the rationale and specifics of simulations in the Lake Sector. Lake model simulations were conducted in two spatial domains: local and global (Fig. 1). Climate change impacts were simulated after model calibration for specific lakes in the local domain (see Sect. 3.1), and for “representative” lakes without calibration in the global domain for each lake-containing grid cell in the ISIMIP global grid (see Sect. 3.2). These two complementary spatial domains balanced the need for site-specific information and the need for a global assessment of climate change impacts on lakes. All temperature simulations were conducted under the assumption that the water level of the lakes remained constant and, therefore, the lakes were decoupled from their watersheds. This assumption allowed us to evaluate lake thermal structure based on meteorological forcing data only and was judged acceptable for the existing phases of ISIMIP (see Sect. 3.5.3). It is anticipated that subsequent simulations, especially those that will evaluate changes in lake biogeochemistry, will abandon this simplification and contribute to cross-sectorial collaborations. On a regional scale, coupled hydrologic and lake model simulations to evaluate changes in lake water level (Hanson et al., 2021) and biogeochemistry (Zwart et al., 2019) for 3692 lakes in northern Wisconsin and Michigan have already been developed. Such regional studies can serve as a model for future cross-sectorial simulation in the ISIMIP.

Figure 1 Map of lakes at local (a) and global (b) scales participating in the ISIMIP2a/b Lake Sector. In panel (a), the local lake sites are visualized through semi-transparent markers; hence darker markers highlight locations where several lakes are located close to each other.​​​​​​​

3.1 Case-study lakes in the local domain

Lakes in the local domain had sufficient information to allow the lake models to be parameterized using individual lake bathymetry and to be calibrated against measured water temperature profiles. Consequently, the local lake dataset was a unique resource for testing and evaluating lake model performance.

For ISIMIP2, bathymetric data and historical data on water temperature from 52 lakes and 10 reservoirs (Fig. 1a, Table S1) were shared among participating modelling teams. Since reservoirs were treated like regular lakes in the simulations, all waterbodies are hereafter called “lakes” (see Sect. 4.3). The geographical distribution of the lakes encompassed a gradient of five major climatic groups in the Köppen–Geiger climate classification, including tropical, arid, temperate, boreal, and polar. Temperate and boreal lakes located in the Northern Hemisphere comprised 87 % of all case-study lakes. The surface area of lakes ranged from 0.011 to 2700 km², with an average and median area of 121.1 and 8.9 km², respectively. Two-thirds of lakes covered surface areas between 1 and 100 km². The average and median mean depths of lakes were 26.3 and 10.8 m (range: 1.7–304.8 m), where 90 % of the lakes were deeper than 3 m. The Secchi depths reported for 49 of the lakes were 4.9 m (average) and 3.5 m (median) and ranged from 0.5 to 32 m, which indicated a wide range of lake trophic status.

In ISIMIP3, the same approach is followed, but the number of lakes has increased. This is achieved through a data call to the research community, and by capitalizing on existing data harmonization efforts (e.g. Pilla et al., 2021).

3.2 Representative lakes in the global domain

Lake simulations in the global domain considered a single, generic lake in each grid cell that contains lakes in the 0.5^∘ × 0.5^∘ ISIMIP global grid. For a given grid cell, such a lake is termed “representative” because it is assumed to represent real lakes bound by its coordinates by sampling their bathymetric information to perform uncalibrated lake model simulations (Subin et al., 2012). The background data and sampling methods for generating representative lakes has evolved from ISIMIP2 to ISIMIP3.

In the global domain of ISIMIP2, generic lakes in each grid cell used average lake depth and surface area information from a rasterized version of the Global Lake and Wetland Database (see Sect. 3.5.1; Lehner and Döll, 2004). Using the mean lake characteristics in the grid cell could not account for the spatial distribution, variability, and non-uniform distribution in the depth and area of lakes in a grid cell. However, similar representations of lakes were used in Earth system and numerical weather prediction models (Subin et al., 2012; Balsamo et al., 2012; Thiery et al., 2015, 2016, 2017; Vanderkelen et al., 2021), and the grid cell lake representation was a necessary trade-off between computational feasibility and global representativeness. This global-scale lake coverage of 17 500 generic lakes made it possible to represent lakes in all major climatic classes and their subclasses, which was not possible in the local lake domain.

In the global simulations of ISIMIP3, this process was improved with a better method of characterizing generic lakes for each grid cell to represent true lakes with morphological characteristics from newer databases: HydroLAKES (Lehner and Messager, 2016​​​​​​​) and Globathy (Khazaei et al., 2022). With this methodology, the 41 449 generic lakes in ISIMIP3 represent true lakes in a more realistic way than for generic lakes defined in ISIMIP2 (see Sect. 3.5.1).

3.3 Lake models participating in ISIMIP

Currently, 10 different lake impact models participate in the ISIMIP Lake Sector, where for some models two different versions were applied (Table 1). There are eight lake models providing calibrated simulations in the local domain: air2water4par, air2water6par, ALBM, FLake, GLM, GOTM5.1, MyLake, and Simstrat. The global ensemble consists of six lake impact models: ALBM, CLM, GOTM5.3, LAKE, Simstrat-UoG, and VIC-LAKE. In the sections below, these models are briefly described. For the models that contribute to both the local and global spatial domains, the global impact model section (Sect. 3.3.2) only describes the differences compared to the local version of the model used.

Table 1 Overview of lake impact models participating in the Lake Sector of ISIMIP2a/b. For “Spatial domain”, L defines local or site-specific and G defines global simulations.

^* Y – Yes; N – No.

3.3.1 Lake models for local simulations

air2water is a hybrid of a physically based and statistical model which simulates lake surface water temperature and epilimnion thickness solely based on air temperature as external forcing (Piccolroaz et al., 2013). The model estimates lake temperature in a single layer characterized by a time-varying thickness according to an empirical relationship accounting for the effect of thermal stratification. Within ISIMIP, two different versions of this model provided simulations for local lakes: air2water4par and air2water6par; these two set-ups of the model differ in the number of parameters affecting the lake thermal dynamics (Piccolroaz et al., 2016). The air2water model has been applied in lakes of varying climatic and morphometric conditions worldwide (Toffolon et al., 2014; Prats and Danis, 2019; Piccolroaz et al., 2020).

The Arctic Lake Biogeochemistry Model (ALBM) is a one-dimensional process-based coupled lake hydrodynamic and biogeochemistry model (Tan et al., 2015). The model simulates water temperature dynamics, ice phenology, and phytoplankton as well as dissolved nitrogen, oxygen, carbon dioxide, and methane. ALBM was originally developed for Arctic lakes (Tan et al., 2015, 2017) but has been used for other lakes across the globe (Guo et al., 2020, 2021; Tan et al., 2018). The thermal regimes of lakes are simulated in ALBM using 1D thermal diffusion equations in both water and sediment columns with atmospheric boundary conditions driven by sensible heat, latent heat, thermal radiation, and solar radiation. ALBM simulates 51 irregular lake layers. Snow and ice dynamics are solved using one snow layer, one white/grey ice layer that is formed when too much snow is accumulated, and multiple black ice layers (Tan et al., 2018).

FLake is a one-dimensional model specifically designed to represent the effects of inland waters in climate models and numerical weather prediction systems (Mironov, 2008). FLake uses a two-layer parametric representation of the lake water column. The upper layer is vertically homogeneous, representing the surface layer produced by wind and convective mixing at the lake surface. The lower layer represents the thermally stratified part of the water column. Two additional layers simulate the ice cover and the lake sediment. The vertical temperature distribution in each layer is modelled by a parameterized function of vertical coordinates, derived from a self-similar representation of the temperature profile. For calculation of surface heat fluxes, the model input includes standard meteorological variables describing the air–lake interaction: air temperature and humidity, wind force, and long-wave atmospheric radiation (or cloud amount for its calculation). The short-wave solar radiation enters the model equations as the volumetric source term distributed across the water column. The FLake version used here uses longwave radiation as a direct input, instead of calculating it from cloud cover.

The General Lake Model (GLM, v3.0) is a one-dimensional hydrodynamic lake model, which simulates temperature stratification in lakes (Hipsey et al., 2019). It uses a flexible Lagrangian grid, and an energy budget approach to simulate mixing. The vertical layer structure can change in number and thickness throughout a simulation, following changes in stratification and lake volume. In this study, we based the initial number of layers on the initial water depth. In addition, GLM includes modules for surface heat exchange and ice–snow dynamics, vertical mixing, and water balance dynamics. GLM can be coupled to the Aquatic Ecodynamics Modelling Library (AED) to simulate water quality dynamics and ecosystem interactions.

The General Ocean Turbulence Model (GOTM v5.3) is a one-dimensional model that simulates the most important hydrodynamic and thermodynamic processes related to vertical mixing (Umlauf and Lemmin, 2005​​​​​​​). GOTM was developed by Burchard et al. (1999) for modelling turbulence in the oceans, but it has been recently adapted for use in lakes (Sachse et al., 2014). 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).

MyLake (v1.12) is a one-dimensional process-based model used to simulate physical, chemical, and biological dynamics in lakes (Saloranta and Andersen, 2007). The model simulates thermal stratification, lake ice and snow cover, and phosphorus–phytoplankton dynamics. It also contains a simple sediment box model. MyLake runs at a daily time step using regularly spaced water layers whose vertical resolution is defined by the user. Different versions of the open-source code have been applied to simulate algal blooms (Moe et al., 2016), CO_2(g) and CH_4(g) (Kiuru et al., 2019), internal phosphorus loads (Markelov et al., 2019), and light attenuation dynamics (Pilla and Couture, 2021).

Simstrat (v2.1.2) is a one-dimensional hydrodynamic model, which specifically includes vertical mixing induced by internal seiches that is not included in most other models (Goudsmit et al., 2002). The model uses layers of fixed depth and supports multiple options for external forcing, comprising several meteorological variables or surface energy fluxes. The model simulates thermal stratification and ice and snow formation (Gaudard et al., 2019). Simstrat has been applied in lakes of varying climatic and morphometric conditions (e.g. Thiery et al., 2014a; Kobler and Schmid, 2019; Mesman et al., 2020) and is operationally applied to provide near-real time, open-access simulation output of the thermal structure and ice cover of all natural Swiss lakes and lake basins greater than 1 km² and a growing number of reservoirs and small lakes (Gaudard et al., 2019).

3.3.2 Lake models for global simulations

The Community Land Model (CLM) Version 4.5 of CLM (Lawrence et al., 2011; Oleson et al., 2013) is a land surface model that includes simulations with the Lake, Ice, Snow and Sediment Simulator (LISSS; Subin et al., 2012). The CLM4.5 model has been used by multiple ISIMIP sectors with one consistent set-up. CLM4.5 simulations and their outputs have been analysed to assess climate change impacts across a range of indicators within ISIMIP (e.g. Schleussner et al., 2018; Lange et al., 2020; Ito et al., 2020; Gudmundsson et al., 2021; Pokhrel et al., 2021; Gädeke et al., 2021; Reinecke et al., 2021; Thiery et al., 2021).

LAKE is an extended one-dimensional model that simulates thermodynamic, hydrodynamic, and biogeochemical processes in the water column and the bottom sediments of the lakes (Stepanenko et al., 2016). The model simulates vertical heat transfer considering the penetration of short-wave radiation (Heiskanen et al., 2015), ice, snow, and bottom sediments. The model explicitly accounts for the exchange of momentum, heat, and dissolved gases between water and the inclined bottom.

VIC-LAKE is a 1D lake model derived from the Variable Infiltration Capacity (VIC) Macroscale Hydrologic Model (Bowling and Lettenmaier, 2010​​​​​​​) and optimized for simulations at a sub-daily timescale. The model is based on a lake energy balance by Hostetler and Bartlein (1990), Hostetler (1991),​​​​​​​ and Patterson and Hamblin (1988). Turbulent mixing is solved with Henderson-Sellers thermal diffusion models using parametrized eddy diffusivity (Henderson-Sellers, 1985). The model also contains an ice module, which dynamically simulates lake ice and ice snow cover.

Simstrat-UoG v1 is based on Simstrat v1.4 and is an earlier version of the model described above. This version uses an earlier snow and ice formulation from Patterson and Hamblin (1988).

3.4 Input data

3.4.1 ISIMIP2a

The Lake Sector simulation protocol for ISIMIP2a was completed in early 2020. This phase focused on the calibration of eight models in the local domain toward projecting these models with meteorological forcings from the gridded ISIMIP2a observations. The same set of gridded meteorological observations was then used for global lake simulations. However, calibrating these lake models globally was unfeasible because of a lack of measured lake water temperatures at this scale.

Meteorological data from 1979 to 2016 at the ISIMIP grid scale (EWEMBI, “EartH2Observe, WFDEI, and ERA-Interim data Merged and Bias-corrected for ISIMIP”, Lange, 2019a) were used as inputs for calibrating the local lake models. From EWEMBI, the grid cell from each local lake's geographical location (Fig. 1 and Table S1) was used for the model calibration. Since the majority of lakes lacked nearby weather stations, the uniform EWEMBI data allowed us to include a broader diversity of lakes and avoid cumbersome data harmonization. Since the EWEMBI dataset was also used to help bias-correct the future climate scenarios used in the ISIMIP2b simulation round (Frieler et al., 2017; Lange, 2019b), the performance of calibrated lake models can be indicative of their ability to simulate past and future climate change when forced by the ISIMIP bias-corrected data. In addition to the calibration runs that were limited to periods when observed water temperature data were available, the local sector modellers were encouraged to drive their lake models with the complete EWEMBI data record between 1979–2016. This was aimed at evaluating the lake models' abilities to reproduce effects of observed meteorologic variability and extreme events on thermal simulations. These simulations could also be used for benchmarking simulations forced with modelled future climate conditions from GCMs. In addition to the EWEMBI dataset, five other reanalysis datasets were provided in ISIMIP2a for modellers to use as inputs according to their capacities, with the goal of exploring the effect of input data choice on simulation outcomes. All datasets are described and referenced in the simulation protocol document (https://www.isimip.org/protocol/2a/, last access: 23 May 2022).

Data providers supplied historical measured water temperature profiles for 62 lakes (Fig. 1a, Table S1). Lake data had to meet two criteria to be included in the local lake dataset: (1) data needed to overlap with the EWEMBI time span, and (2) temperature profiles needed to encompass at least two consecutive years in the case of (sub)-daily sampling frequency, or at least five consecutive years in the case of (sub)-monthly sampling frequency. These criteria enabled intra- and inter-annual variability to be captured in water and meteorologic conditions in the model calibration procedure. A few lakes from under-represented geographical locations (e.g. tropics) were included despite shorter water temperature records.

Water temperature data were harmonized to a uniform data format and visually quality-controlled to remove outliers. In addition to water temperature, the data providers supplied detailed information of the lake depth and hypsometric data to characterize the morphometry of lake basins, which are required as input to most of the lake models.

3.4.2 ISIMIP2b

ISIMIP2b was designed to compare lake responses to simulated historical and projected future climates relative to pre-industrial climates with a focus on improving the understanding of the effects of global warming in the range of 1.5 to 2 ^∘C (Frieler et al., 2017). The lake ensemble included simulations forced by bias-adjusted data from four GCMs, covering historical and up to three representative greenhouse gas concentration pathways (RCPs): a low (RCP2.6), medium–high (RCP6.0), and high emission scenario (RCP8.5). The past and future responses of lakes from these simulations were compared to simulations forced by bias-adjusted pre-industrial control (picontrol) climate data from the same four GCMs to quantify differences from pre-industrial conditions. These differences can be thought to represent the “pure” effect of ongoing changes in climate on simulated lake water temperatures, with minimal confounding effects from changes in further human influences that were identical between the two sets of simulations.

Climate input data for ISIMIP2b were derived from four GCMs, namely GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, and MIROC5 (Taylor et al., 2012; Frieler et al., 2017) that were available from CMIP5. These GCMs were chosen since they best met the needs of all sectors participating in the ISIMIP and provided the necessary scenario length at daily temporal resolution. They had a wide range of projected warming rates, with GFDL-ESM2M and HadGEM2-ES representing the lower and higher ends of the warming spectrum, respectively. The data management team at the Potsdam Institute for Climate Impact Research (PIK) bias-adjusted the GCM data with a reference dataset of atmospherics observations (EWEMBI, Lange, 2019a) using the statistical transfer functions by Hempel et al. (2013) modified to correct known biases in modelled variables (Frieler et al., 2017). The bias correction was aimed at preserving trends and distributions of modelled variables relative to observed atmospheric observations. All meteorologic data except horizontal wind components were bias corrected. The list of output meteorological variables from GCMs that were used to drive the lake models can be found in Table 2.

Table 2 Standard climate forcing variables used to govern lake models to simulate past, present, and future changes in local and global lakes. Dots between brackets indicate optional input variables.

^* Uncorrected.

Lake Sector simulations followed the ISIMIP2b protocol (Fig. 1, Frieler et al., 2017, https://www.isimip.org/protocol/#isimip2b, last access: 23 May 2022). To estimate the effects of historical climate warming, lake model simulations forced with data from historical emission scenarios were compared with simulations forced with data from the picontrol scenario. Likewise, to evaluate future climate impacts, lake model simulations were forced with data from the RCP trajectories (RCP2.6, RCP6.0, and RCP8.5) and compared to results from simulations forced with picontrol data. The time spans of different climate change scenarios were 1661–1860 (picontrol), 1861–2005 (picontrol and historical), and 2006–2099 (picontrol, RCP2.6, RCP6.0, and RCP8.5). An extended period between 2100 and 2299 was also used for simulations based on available results for specific emission scenarios (picontrol, RCP2.6) to evaluate longer-term changes in global temperature that meet the Paris Agreement objectives (Frieler et al., 2017).

Since the lake modelling strategy was specifically designed to evaluate only thermal changes resulting from changing atmospheric conditions under the assumption of no watershed inputs (constant lake level), the simulations are not influenced by any changes in land use or socio-economic conditions that would affect watershed inputs to the lakes or changes driven by changes in lake trophic status. The pre-industrial reference simulations (picontrol) assumed fixed socio-economic conditions and land use (1660–1860). CLM4.5 provides additional sets of simulations according to protocols for other ISIMIP sectors (i.e. biomes, agriculture, water (global), permafrost) with a combination of socio-economic (1860-soc and 2005-soc) and CO₂ fertilization (2005-co2) scenarios. Lake temperature simulations that were a component of CLM4.5 did not account for these additional scenarios.

3.4.3 ISIMIP3a/b

The protocol for the third simulation round of ISIMIP, ISIMIP3, is currently ongoing. ISIMIP3 is largely similar to the ISIMIP2 protocol but includes counterfactual climate forcing in ISIMIP3a, and the next generation (CMIP6) of climate model forcing and various emission scenarios in ISIMIP3b. Below, we highlight the main differences between both simulation rounds.

In ISIMIP3a, the observational climate forcing covers the period 1901–2016 and consists of the Global Soil and Wetness Project version 3 (GSWP3; Dirmeyer et al., 2006), homogenized to W5E5 for the period 1901–1978 (Lange, 2019a) and a combination of the W5E5 dataset (Lange, 2019b; Cucchi et al., 2020) for the period 1979–2016 and GSWP3 before that. This observational dataset was bias-corrected using observations from the Global Precipitation Climatology Centre (GPCC) (more details at https://climatedataguide.ucar.edu/climate-data/gpcc-global-precipitation-climatology-centre, last access: 23 May 2022) and Climatic Research Unit (CRU) (more details at https://crudata.uea.ac.uk/cru/data/hrg/, last access: 23 May 2022) using the method outlined in Lange (2019b). In addition to providing data for the calibration for local lake models, ISIMIP3a provides counterfactual climate forcing, which is a detrended version of the historical climate forcing (Mengel et al., 2021). Models driven by the counterfactual climate and other historical human pressures provide a baseline to compare with simulations forced by the observational climate forcing to determine climate change impacts, paving the way for IPCC Working Group II style impact attribution (Cramer et al., 2014). In ISIMIP3b, the climate forcing is updated to include the next generation of CMIP6 simulations, which were bias-corrected with a new adjustment routine correcting the simulations towards the W5E5 observational data (https://www.isimip.org/gettingstarted/isimip3b-bias-correction/, last access: 23 May 2022). Climate simulations from five GCMs were provided, namely GDFL-ESM4, UKESM1-0-LL, MPI-ESM1-2-HR, IPSL-CM6A-LR, and MRI-ESM2-0. In addition to the picontrol and historical emissions simulations, future simulations include the SSP1-RCP2.6, SSP3-RCP7, and SSP5-RCP8.5 emission pathway scenarios. Following the CMIP6 protocol, the simulation periods were updated to 1601–1849 for pre-industrial, 1850–2014 for historical, and 2015–2100 for future simulations. Like in ISIMIP2, the GCMs chosen for ISIMIP3 were constrained by data availability, yet they are also a subset of better-performing models relative to the entire CMIP6 ensemble and they contain structurally independent model components (Lange, 2021).

3.4.4 Climate data application

All bias-adjusted meteorologic forcing data provided by ISIMIP sectors have a daily temporal resolution and a spatial resolution matching the ISIMIP grid scale. While most models in the Lake Sector performed simulations at daily time steps, some models required temporally disaggregated forcing data at sub-daily time steps. The modelling teams performed temporal disaggregation using their customary approaches (see Sect. 4). For simulations in the local domain, data were extracted for grid cells corresponding with the lakes' geographic locations. No further downscaling or local corrections were applied to ensure consistency in the forcing data applied to all local lakes.

3.5 Lake parameterization

In ISIMIP2 and ISIMIP3, to account for variations in individual lake responses to meteorological drivers (Heiskanen et al., 2015; Kraemer et al., 2015; Shatwell et al., 2019), there were only two types of data needed by the lake models: a description of the lake bathymetry and information on the lake water transparency, which are necessary for estimating the diffuse attenuation coefficient of incoming shortwave radiation.

3.5.1 Bathymetry

Most lake hydrodynamic models require the hypsographic relationship between depth and surface area, which is critical for determining layer volumes and storage and the vertical transfer of heat. Data providers supplied these bathymetric data for each lake in the local domain. The two versions of air2water models did not require information on lake bathymetry.

For global lake simulations in ISIMIP2, the bathymetry of the representative lakes in each grid cell was derived from a rasterized version of the Global Lake and Wetland Database (GLWD; Lehner and Döll, 2004; Toptunova, 2019). Specifically, for each grid cell, average lake depth and lake surface area values were calculated from all GLWD lake data contained within the grid cell. Lake bathymetry for each generic lake then was assumed to be cylindrical. In the case of the LAKE model, both surface area and mean depth of global lakes were obtained from GLDBv2 (Choulga et al., 2014). In the case of the global simulation model CLM4.5 (see Sect. 3.3.2), all representative lakes had a constant 50 m depth. For all of the lake models, a cylindrical shape was assumed to represent lake bathymetry. The gridded lake masks for the surface area and mean depths can be accessed online (https://data.isimip.org/, last access: 23 May 2022).

In the ISIMIP3 global lake simulations, we selected a representative lake for each grid cell from the 1.4 million lakes included in the HydroLAKES shapefiles (Lehner and Messager, 2016​​​​​​​). Here, each lakes was assigned to a single grid cell using the location of the lake centroid. For grid cells containing more than one lake centroid, we selected the lake with the depth corresponding to the depth weighted median (weighted by the area of the lakes) of all the lakes contained in the pixel. Then, for each representative lake selected in the previous step, volume, area, mean and maximum depth, and 11-level area and volume hydrographic curves were extracted from GLOBathy (Khazaei et al., 2022). After applying this procedure, we obtained a total of 41 449 representative lakes across ISIMIP grid cells (see Fig. 3).

3.5.2 Water transparency

Water transparency may mediate a lake's response to climate change (Butcher et al., 2015, Magee et al., 2016; Shatwell et al., 2019). The attenuation coefficient (K_d, m⁻¹) of shortwave radiation is a key parameter to describe water transparency in the lake models (Potes et al., 2012). In all simulations, single non-varying values of K_d were used.

For simulations in the local domain, multi-season and multi-year water transparency data were available from 49 lakes that were used to calculate a diffuse attenuation coefficient (K_d) value or a mean Secchi depth (Z_SD, m) (Table S1). When only Z_SD measurements were available, K_d was estimated at $\mathrm{1.7}/Z_{\mathrm{SD}}$. (Poole and Atkins, 1929). If both mean Z_SD and K_d were provided, the directly measured K_d was used.

When both K_d and Z_SD were lacking, we approximated Kd as a function of mean lake depth following the equation derived for 88 Swedish lakes (Håkanson, 1995​​​​​​​),

$$K_{\mathrm{d}}=\mathrm{1.1925}\times max(\text{mean\_depth},\mathrm{1}{)}^{-\mathrm{0.424}},$$

or from maximum depth following the expression derived for 1258 global lakes (Woolway et al., 2021b​​​​​​​):

$$K_{\mathrm{d}}=\mathrm{5.681}\times max(\text{max\_depth},\mathrm{1}{)}^{-\mathrm{0.795}}.$$

For ISIMIP2a, the different modelling teams defined the best method to parameterize transparency based on the specific lake model requirements and previous protocols developed for calibration and simulation with any given model. Consequently, transparency parameterizations varied both with lake models and the spatial domain of simulations. In the local domain, K_d values derived from Z_SD or Hakanson's expression were used in ALBM, Simstrat, and GLM, whereas FLake runs adopted the approach outlined in Woolway et al. (2021b). In GOTM and MyLake, the mean K_d was determined in the calibration process. The two air2water models did not require water transparency. In the global-scale simulations, the grid-varying mean depths of the lakes were used to estimate K_d values from Hakanson's expression in all lake models except for VIC-LAKE. This estimation process is also in the CLM4.5 (Oleson et al., 2013). In the VIC-LAKE model, two-band (visible and near-infrared) Beer's law radiation constants were used to parameterize transparency (Bowling and Lettenmaier, 2010).

3.5.3 Water balance

To simplify lake simulations, the water balance and water inputs and withdrawals were not considered in ISIMIP2 and ISIMIP3. The formulations of some lake models (e.g. air2water or FLake) do not explicitly include hydrological balances. For the rest of the models, the precipitation and evaporation component of water mass exchange was switched off (i.e. only heat exchange occurred) or compensated with a closure term (e.g. CLM4.5). This assumption allowed us to evaluate changes in lake thermal structure in the time frame of the ISIMIP2 and ISIMIP3 simulation periods.

Regional studies assessing the hydrologic responses of lakes to an ensemble of future climate change scenarios show that our omission might variably affect lakes depending on lake type and future climate outcomes for seasonal drying and wetting (Hanson et al., 2021; Hunt et al., 2013). These studies found that drainage lakes in northern Wisconsin, United States, which are hydrologically mediated by lake inflows and outflows, were projected to maintain stable water levels because of competing climatological factors that did not promote a clear drying trend. Under our omission of lake water balances, projections for such lakes could lose reliability where future climate conditions reduce watershed runoff. In the same region, seepage lakes with minimal surface water fluxes and a greater dependence on groundwater inflows, however, were projected to significantly decrease in water level, especially in higher-elevation regions near groundwater divides. These studies are relevant for both our local and global lake simulations. For lakes in the local domain, despite accurate representations of historical changes in lake thermal structure (Table 3), the omission of a water balance could additionally affect the simulated climate change impacts in seven lakes and reservoirs with large water level fluctuations (Table S1); thus caution should be used when evaluating these results. For lake simulations in the global domain, this omission is yet another necessary trade-off between experimental complexity and spatial representativeness (see Sect. 3.2).

Table 3 Summary of local domain calibration of lake models participating in the ISIMIP2a/b simulation rounds. Summary statistics for calibrated parameters and models' performance represent the cross-site mean, median [minimum – maximum].​​​​​​​

^a a_seiche parameter was calibrated for lakes <100 km⁻² (n=51). ^b a_seiche_s and a_seiche_w parameters were calibrated separately for summer and winter months in lakes >100 km⁻² (n=8).

3.6 Calibration of local lake models in ISIMIP2a

Eight lake models had specific parameters and coefficients calibrated based on what each modelling group felt was appropriate for use with their specific lake model (Table 3). Each modelling group defined reasonable coefficient ranges based on past experience and the physical constraints that would set limits on the parameter and coefficient values. For each model, the same calibration routine and objective function was applied to all lakes in the local domain. Different objective functions (e.g. RMSE, NSE, Pearson r​​​​​​​; see Table S2) were adopted by the different models so that modellers could use their optimal criteria for calibration. In all cases the model performance was optimized by minimization of the difference between simulated and measured water temperature.

The number of calibrated parameters and coefficients in a specific model ranged from one (FLake) to nine (ALBM, Table 3). The calibrated coefficients were mostly related to processes controlling surface heat and energy fluxes, turbulent kinetic energy and wind stress, and light attenuation. Other calibrated coefficients for specific processes were model-specific, including sediment structure and heat fluxes (ALBM), seiches (Simstrat), and ice–snow energy fluxes (MyLake, Simstrat). To allow for comparing the lake models' performance in predicting measured water temperature, for all lake models two common metrics of model fit were calculated in post-processing (not necessarily coincident with the calibration metrics): the root-mean-square error (RMSE) and coefficient of determination (R², Table 3).

Most lake models were calibrated with the full series of available measured observations. In this majority of cases, no data were withheld for an independent model validation. Considering the relatively short temporal extent of lake measurements, this was done to base parameter estimates on the full range of environmental conditions encountered during simulation for producing robust future projections (Larssen et al., 2007). This is justifiable given extensive research validating the performance of these models outside the calibration period (e.g. Stepanenko et al., 2013; Thiery et al., 2014a) and arguments calling for scepticism of the split-sample approach to calibration and validation (Augusiak et al., 2014; Shen et al., 2022). Exceptionally, ALBM only used the full series of measured observations when the observations were shorter than 5 years. Where measurements exceeded 5 years, modellers running ALBM simulations opted for a split-sample approach to tuning their model and used the first 5 years of measurements for calibration.

3.7 Long-term simulations in ISIMIP2b

The ensemble of lake models in both the local and global domains was forced with the bias-adjusted GCM outputs for the pre-industrial control, historical, and future climate change scenarios. When running the long-term simulations, the calibrated models for each local lake were used, so that each model was optimized for that lake based on the historical calibration described in Sect. 3.6. Spin-up periods used with the local lake models varied and were dependent on the protocols and experience of each modelling group (Table S2). When a spin-up period was used, the spin-up data were created either by repeating the initial year(s) of the scenario input data and then adding these duplicate data to the beginning of the forcing data, or by using a portion of the historical scenario to spin-up future scenario simulations. Initial conditions used for water temperature profiles in the local lake simulations also varied with model and geographical location and were based on either observed temperature profiles, an assumed isothermal 4 ^∘C profile, or related to the mean annual air temperature at the local lake location. A more detailed description of the modelling workflows that were used to spin up and initialize each model in the local domain is given in Table S2.

For simulations in the global domain, most lake models used default parameter and coefficient values that were set according to previous experience with each model (see Table 1: “Key references”). Exceptionally, for GOTM, the average values of calibrated coefficients from the GOTM local lakes (Table 3) and default values for the coefficients that were not calibrated (Umlauf and Lemmin, 2005; Sachse et al., 2014) were used for all representative lakes in the global domain. Similarly, the methods for spinning up the models in the global simulations also varied depending on the practices applied by each modelling group. Groups working in the global domain tended to use longer spin-up periods and data from either the picontrol or historical scenarios to create the spin-up data that were added to the scenario forcing data. The initial water temperature profiles used in the global lake simulations also varied. In some cases, the models were initialized as homogenous profiles often based on the mean annual air temperature or linear profiles based on the mean annual air temperature and an assumed 4 ^∘C bottom temperature, or linear profiles using a fixed surface and bottom temperature. A description of the modelling workflows that were used to spin-up and initialize each model in the global domain is given in Table S3. More detailed model-specific simulation set-ups can also be found at https://www.isimip.org/impactmodels/ (last access: 23 May 2022).

3.8 Output data format

All outputs from the models were aggregated to daily averages (Tables 4, S4). The vertical resolution of the simulated water temperature profiles in the local domain was reported at 0.5 m intervals for lakes with <50 m maximum depth and at 1 m intervals for lakes >50 m. However, the vertical resolution of simulated temperatures in the global domain was limited by file storage capacity. The number of reported layers depended on the depth of the representative lake and ranged from 1 to 13 (GOTM, LAKE, Simstrat-UoG), from 1 to 50 (ALBM), and from 1 to 1000 (VIC-LAKE). Output from CLM4.5 was grid-invariable, representing water temperature in 10 layers. The remainder of reported variables (thermodepth, surftemp, bottemp) represented a single value, which was either calculated using the approach presented in the simulation protocol (see https://www.isimip.org/protocol/2b/, last access: 23 May 2022)​​​​​​​ or was directly outputted by the lake model (Tables 4, S4). The Lake Sector simulation protocol provides the model performance metrics used during calibration of lakes in the local spatial domain (Table 3). A full list of variables simulated within ISIMIP2b is summarized in Table S4.

Table 4 Common output variables reported by local (L) and global (G) models participating in the Lake Sector of ISIMIP2a/b. The variable watertemp is a full water temperature profile. Naming of lake models and variables are ordered in an alphabetical order (see Table S4 for a list of full variable names).

^* The model provides a time-varying estimate of the well-mixed surface layer participating to the heat exchanges with the atmosphere.

This diversity of GCM input datasets, emissions scenarios, lake models, and their output variables means that the total ensemble of impact simulations under the Lake Sector requires considerable storage space. In addition, those appreciable computing resources should be anticipated by potential future collaborators. For example, the global lake simulations for ISIMIP2b take up 14 TB of storage space. This means that applications with simulations under multiple GCMs, lake models, and scenarios for a given variable will require high-performance computing resources. For running simulations, computing times may vary depending on the scale of one's contribution. On the one hand, simulating a local, calibrated lake with FLake for a single scenario and GCM combination may take seconds on a laptop (Thiery et al., 2014b), but, on the other hand, global simulations from CLM4.5 for one such scenario and GCM combination will require several weeks using 144 compute cores on a high-performance computer, substantiating both computational costs and resources for dataset storage. These technical prerequisites, in addition to individual model feasibility issues for local versus global domain simulations, explain the discrepancy in model availability across the ISIMIP2 local and global simulations.

4 Results and discussion

4.1 Results from ISIMIP2a

Calibration and performance of local lake models

The​​​​​​​ simulated water temperatures from the calibrated lake models compared well to the measured water temperature data using the entire record of recorded water temperature data from each local lake. Based on the simulation data from 62 lakes, all eight local lake models were calibrated with a multi-model mean RMSE of 1.50 ^∘C that ranged from 0.98 (air2water6par) to 2.41 ^∘C (FLake, Table 3). The coefficients of determination (R²​​​​​​​) ranged from 0.59 (MyLake) to 0.96 (air2water6par), with the multi-model average value of 0.84. The lake models predicting surface water temperature only, air2water4par and air2water6par, showed lower prediction errors compared to the lake models predicting full water temperature profiles. While the multi-model mean goodness of fit was reasonable for most lakes, 16 % of lakes showed RMSE larger than 2 ^∘C, indicating less certain predictions (Tables 3, S1). For individual models, the number of lakes with RMSE exceeding 2 ^∘C varied from 3 lakes (5 %, air2water4par) to 40 lakes (65 %, FLake).

Although the ISIMIP2a forcing data used a daily time step at the ISIMIP grid scale, the prediction errors in water temperature were relatively small (Table 3), even though these input data and their resolutions are, in general, less than optimal for the simulations of individual lakes (Bruce et al., 2018). An exception is the air2water model that, owing to its statistical and data-driven calibration of its model's parameters, has been shown to be able to provide the same projections irrespective of the nature of the air temperature dataset used to drive the model (Piccolroaz et al., 2018). It should be noted, however, that inter-model performance comparisons are difficult here. Due to the diverse discretization of lake temperature profiles across models, each model is being evaluated on a derivation of available lake measurements. Therefore, the observations used as a reference in the performance metrics are different across models, reducing the comparability of model performances for a given metric. The average errors in the prediction of water temperature observations were comparable with previous multi-model and/or multi-site modelling studies, where the mean RMSE in water temperature predictions ranged from 1.10 to 2.79 ^∘C (Stepanenko et al., 2013; Winslow et al., 2017; Bruce et al., 2018; Piccolroaz et al., 2020). Similarly, the prediction of epilimnetic temperature showed lower errors compared to predictions of hypolimnetic temperature (Winslow et al., 2017; Bruce et al., 2018).

4.2 Results from ISIMIP2b

Impacts of past and future climate change on lakes

Time series​​​​​​​ of ensemble simulations of lake surface temperature for local lakes over the historical (1851–2005) and future (2006–2100) periods are shown in Fig. 2a. Each ensemble combines the results from 62 well-studied lakes and three separately calibrated lake models. Mean annual surface water temperatures increased by 0.15 ^∘C at the end of the historical period (present-day, 1976–2005) relative to the pre-industrial control. These simulations support in situ observations showing that lakes across the globe are already warming (Woolway et al., 2020; O'Reilly et al., 2015). Future projections (2006–2099) accounting for low (RCP2.6) to high (RCP8.5) greenhouse gas emissions under present-day socio-economic conditions, provide ensemble estimates of lake surface warming of 1.38, 2.46, and 3.85 ^∘C by the end of the century (2070–2099) relative to pre-industrial control, respectively (Fig. 2a). For example, the strong mitigation measures associated with RCP2.6 resulted in lake surface temperature remaining below 2 ^∘C. The ensemble projections consistently show a slower surface water warming under RCP2.6, starting mid-century (Fig. 2a), than for other greenhouse gas emission scenarios.

Figure 2 Example output data in the ISIMIP2b Lake Sector. (a) Local model time series of mean annual surface temperature derived from ensemble simulations for a period 1861–2100 (historical and future) averaged across lakes and climate models. Colour-coding indicates greenhouse gases scenarios; line types differentiate lake models. Warmer temperatures simulated by the GOTM model are due to the fact that the model version used had only a very rudimentary method for simulating the effects of lake ice on the surface heat exchange. Spring–fall GOTM temperature simulations were compatible with the other models. (b) Temperature changes by 2070–2099 compared to pre-industrial levels in summer mean temperature at two lake depths simulated with GOTM (local) model evaluated with data from four climate models under RCP8.5. (c) Global outlook for the mean annual surface temperature of lakes by the end of the 21st century (2070–2099) under RCP2.6 simulated with GOTM global.

Figure 3 Maps of lake area (a), lake volume (b), and mean lake depth (c) used in the ISIMIP3 simulation round. In the input data for ISIMIP3 simulations, a single lake is assigned to each grid cell. However, here we show a modified version of the dataset to delineate large lakes in the global map. The maps are derived from HydroLAKES (Lehner and Messager, 2016) and GLOBathy (Khazaei et al., 2022) datasets using the ISIMIP3 lake mapping methods described in “Code and data availability”.

Based on the anomalies between the pre-industrial control and future scenarios, all lake models showed similar warming rates and trajectories of change. However, the GOTM local simulations were on average 1.75 ^∘C warmer than simulations from ALBM and Simstrat, probably because the version of GOTM used for the local lake simulations had only a very rudimentary description of the effects of lake ice on surface heat exchange. These results show the importance of using an ensemble of models to increase the robustness of simulated past and future changes and making interpretations less dependent on a single or small suite of the lake models used (Trolle et al., 2014).

The common fundamental output from all lake models was water temperature, and for most models, this output is in the form of a full vertical profile at a daily time step. These data and other related model output were also aggregated to different metrics describing lake hydrodynamics, e.g. thermocline depth, the onset of stratification, energy and heat fluxes at the air–water interface, and lake ice characteristics and dynamics (Table 4). The methods to calculate these metrics are defined in the Lake Sector protocol (https://www.isimip.org/protocol/, last access: 23 May 2022), and for additional metrics, the full lake water temperature profiles can be further processed by users.

The average changes in surface and bottom water temperature for the 62 lakes for the RCP8.5 greenhouse gas emission scenario, using the GOTM model forced by the four GCM outputs, are shown in Fig. 2b. Results from the ensemble simulations of the local lakes' future responses show faster warming of surface waters (local-lake mean 4.08 ^∘C) than bottom waters (1.49 ^∘C) by 2070–2099. On average, the difference between the surface and bottom water temperature anomaly was 2.6 ^∘C. There was a wide range of lake responses in the local domain (Fig. 2b), with an average range in the change in surface temperatures anomalies derived from the ensemble of 2.28 ^∘C and bottom temperature of 3.22 ^∘C. These results are consistent with previous findings of the diverse responses in lake surface temperature across the globe (O'Reilly et al., 2015; Pilla et al., 2020) depending on a complex interaction of climate regions (Piccolroaz et al., 2020), lake morphology (Toffolon et al., 2014), and atmospheric conditions (Spence et al., 2013), and changes in the responses of bottom temperatures being influenced by the lake's morphometric characteristics (Kraemer et al., 2015). The Lake Sector local domain provides information on the lake-specific characteristics related to morphometry and water transparency (Table S1) to enable investigation of how the observed differences in responses to climate warming are influenced by lake characteristics.

The variability in Fig. 2b is a result of both variable lake responses and the differences in the forcing associated with the four GCMs. The mean change in surface water temperatures under RCP8.5 until the end of the 21st century ranged from 2.39 ^∘C (when forced by GFDL-ESM2M) to 5.34 ^∘C (when forced by IPSL-CM5A-LR). A similar pattern was observed for bottom temperature, although the differences were less pronounced (1.19–1.78 ^∘C). The changes in the mean surface temperature followed the differences in the air temperature projected by the four GCMs. Sorted from colder to warmer based on simulated impacts on air and water temperature on the local lakes, the GCMs are ranked in the order of GFDL-ESM2M, MIROC5, HADGEM2-ES, and IPSL-CM5A-LR. Similar differences were observed in water temperature and ice changes by Woolway and Merchant (2019). These results indicate that the choice of the GCM has a large effect on the changes predicted by the lake models. Using outputs from several GCMs, following the ISIMIP protocols, therefore, provides the advantage of including ensemble forcing data in simulations of climate change impacts on lakes, increasing the robustness of predictions (Grant et al., 2021).

The results of global domain simulations made with the GOTM model are shown in Fig. S1 for three greenhouse gas emission scenarios and as an ensemble of four GCMs. Under RCP2.6, the emission scenario with the strongest mitigation measures, the global mean annual lake surface temperature was projected to be 12.7 ^∘C (range: 3.8–29.4 ^∘C) by the end of the 21st century (Figs. 2c, S1a). However, global mean lake surface temperatures of 13.4 ^∘C (4.2–30.4 ^∘C) and 14.3 ^∘C (4.8–31.6 ^∘C) were projected for the medium–high emission (RCP6.0) or high-end emission scenario (RCP8.5), respectively (Fig. S1c, e). Mean annual lake temperature was projected to increase by 0.9 ^∘C (0.53–1.32 ^∘C), 1.7 ^∘C (1.0–2.3 ^∘C), and 2.6 ^∘C (1.6–3.6 ^∘C) under these three greenhouse gas emission scenarios relative to the pre-industrial control (Fig. S1b, d, f). Simulations in the global domain allowed the documentation and visualization of spatial variations in lake thermal structure that are not possible using geographically constrained data in the local domain. The most pronounced spatial pattern was a latitudinal gradient of warming of global lakes (Figs. 2c, S1). These results corroborate previous global-scale modelling studies, although here the results are based on ensemble simulations compared to single model simulations.

Existing studies applying these simulations demonstrate the many possibilities for exploring the impacts of climate change on lake physics under the ISIMIP protocol. ISIMIP simulations have been used in a first ever assessment of the global heat uptake by inland waters (Vanderkelen et al., 2020), a relevant addition to existing evaluations of Earth's global heat budget in its land, atmosphere, and oceans. The ISIMIP Lake Sector database has also been used to assess present and future alterations of lake mixing regimes (Woolway and Merchant, 2019) and the shifts in lake stratification and their climatic drivers (Woolway et al., 2021b). Finally, both event and trend attribution of lake heatwaves (Woolway et al., 2021b​​​​​​​, 2022) and lake ice cover changes (Grant et al., 2021), respectively, have been undertaken using ISIMIP simulations in combination with global observational datasets to confirm the role of anthropogenic climate change in observed lake changes.

4.3 Future work

4.3.1 Model response to observational vs. simulated forcing data

In addition to simulations using ISIMIP2a forcing, the ALBM and FLake models were also used for simulations forced by EWEMBI observational data (1979–2016). This will allow for assessment of the difference in model output when used with observational forcing data compared to simulations with GCM forcings during the historical time period. Given that impacts under past and future climates are modelled with bias-adjusted GCMs, a comparison with simulations using the observed data used for bias correction will allow an assessment of how simulations forced with the GCM historical inputs compare with those forced using observed (historical) climate (see also Piccolroaz et al., 2018 for a similar analysis). This can give an estimate of the uncertainty in the ISIMIP GCM scenarios and the bias correction method. There are so far no studies for this application of the ISIMIP2a simulations, but the existing simulation output archives are publicly accessible and hold potential for further study.

4.3.2 Lake hydrology and water quality assessments

Current lake modelling activities in ISIMIP are biased towards lake physics and concentrate mostly on water temperature and related variables like ice cover or stratification state. Lake managers require more than that and are usually highly interested in projections for water quantity (i.e. inflow discharges) and water quality and potential effects on the services lakes deliver. Future directions of the lake sector beyond ISIMIP3 are therefore seen in (i) linking the water sector (hydrological models) with the lake sector in order to integrate water quantity projections into lake simulations and (ii) adding water quality descriptors by biogeochemical modelling of lakes. Such modelling can make projections for the future development of ecosystem services and biodiversity in lakes in relation to climate change and socio-economic development. For climate change, such assessments can be directly built on the ISIMIP2 and ISIMIP3 simulations of the Lake Sector but require linkage with the transport of water and nutrients from their catchments (Janssen et al., 2019). For that, nutrient transport models such as IMAGE-GNM (Beusen et al., 2015) or MARINA (Strokal et al., 2016) need to be aligned akin to the ISIMIP approach.

4.3.3 Global-scale calibration and validation

Important steps can be made in the development of a global-scale dataset for calibration and validation purposes. There are a few challenges to overcome in the future (Janssen et al., 2015). First, due to project-based research, long-term measurements are rare as often measurement campaigns stop when projects are over. Second, data are often locked within institutes, meaning that a consistent global database requires corporation between various parties. Similarly, in situ data that have not been properly indexed and stored, sometimes referred to as “dark data”, require rescue efforts to extend back our measurement period of lakes. Third, data are gathered inconsistently, for instance by using different methods, measuring over different periods, or collecting at different spatial scales. Remote sensing could overcome these issues to some extent, as they can provide long-term, global observations.

Remote sensing datasets for lakes are increasing (Dörnhöfer and Oppelt, 2016). Examples of already existing datasets are datasets for lake temperature (Sharma et al., 2015), ice phenology (Wang et al., 2021), and even biological indicators (Fang et al., 2022; Hou et al., 2022). A disadvantage is that remote sensing is limited to proxy values, which still require ground truthing by in situ monitoring data. Moreover, remote sensing performs variably depending on the measurement system and weather conditions and are variable in assessment. While optical imagery is easily obscured by cloud cover, active microwave systems can be used in all-weather conditions for some variables such as ice cover (Kilic et al., 2018; Murfitt and Duguay, 2021). Therefore, satellite observations must be combined with highly spatiotemporally resolved in situ measurements from buoys, field sampling programs, and long-term monitoring networks (Rand et al., 2022). Specifically, in situ measurements are essential for observing lake processes below the water surface (such as stratification and mixing), to improve understanding of complex air–water energy fluxes (such as evaporation), and to maintain long-term perspectives that began prior to the advent of satellites and regardless of weather conditions that adversely impact some satellite measurements. First attempts at such databases are for example the HydroLAKES database, which already has water discharge into lakes (Messager et al., 2016​​​​​​​).

5 Conclusions

Modelling the impact of climate change at a global scale using an ensemble of lake models requires data, vision, ambition, and a strong collaborative network of researchers from a range of disciplines. The first ensemble model simulations in the Lake Sector of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) has leveraged such a network to design and execute a protocol that has now provided state-of-the-art scientific evidence of climate change impacts on lakes under low to high greenhouse gas emission scenarios. The Lake Sector protocol in the local domain allows for the calibration of lake models forced with historical ISIMIP2a inputs and parametrized using site-specific data of bathymetry. Comparison of simulated and measured water temperature from 62 well-studied lakes facilitated detailed calibration and evaluation of the models in the local domain. In future global simulations, these locally derived parameter and coefficient values could improve the full ensemble of models that have so far been uncalibrated in their global domain applications. Simulations in the global domain provided daily outputs of lake thermal conditions at the ISIMIP grid scale.

Our simulations at both the local and global resolution quantify past and future changes in water temperature, energy and heat fluxes, and ice with unprecedented geographical coverage. Simulations by the ISIMIP Lake Sector also provide previously unattainable opportunities to evaluate the levels of uncertainty in simulations related to the differences in forcing data between reanalyses, GCMs, emission scenarios, and in model structure and parameterization among lake models.

Simulations by the ISIMIP Lake Sector continue to estimate the range and robustness of plausible lake responses to global warming either at 1.5 ^∘C above pre-industrial levels as defined by the ISIMIP2b protocol or for any other future greenhouse gas emission scenarios. This work furthers the state of the art in freshwater science (Vanderkelen et al., 2020; Grant et al., 2021; Woolway et al., 2021b).

Here, we have described the protocol of the Lake Sector in ISIMIP2 and ISIMIP3, which includes the simplifying assumption that hydrologic inputs from the lake watershed had minimal effects on the simulated thermal structure. While this is a reasonable assumption for lake hydrodynamic simulations, it will not be sufficient for simulations of lake biogeochemistry and ecology that strongly depend on the nutrient inputs from the lake watershed. Under the ISIMIP framework's provision of consistent climate forcing datasets and scenarios, the climate change impacts simulated in the Lake Sector are comparable with simulation results from other ISIMIP sectors, supporting cross-sectoral assessments of climate change impacts (Lange et al., 2020; Vanderkelen et al., 2020; Thiery et al., 2021). Ultimately, we expect that the improved simulations of lake hydrodynamics presented here will form a basis for more complex simulations of water quality, lake level fluctuations, and other greenhouse gas emissions scenarios in upcoming simulation rounds, where lake water quality models can be coupled to the hydrologic and biogeochemical outputs from other sectors of the ISIMIP.