Comparison between methods to predict climate change impacts on tropical shallow lakes

Climate change affects a lake's hydrodynamics, creating consequences for whole ecosystem processes (Amorim et al. 2019; Lewis et al. 2019; Woolway et al. 2019; Piccioni et al. 2021; Sisay Kebede Balcha et al. 2023). Lakes' and reservoirs' behaviour is influenced by morphometric, meteorological, and hydrological conditions, which are directly reflected in the reservoir water quality. Hydrodynamic transport describes how chemical or biological material is moved from one location to another, along with salinity, temperature, and sunlight influencing the kinetic processes in the aquatic environment (Imberger 1998; Amorim et al. 2019; Polli & Bleninger 2019; Merino-Ibarra et al. 2021).

Natural and artificial lakes are an important part of the aquatic ecosystem. They have many uses, like water supply, hydropower generation, irrigation, flood control, and a landscape component. In many cases, they have multiple uses. Additionally, specific hydrodynamic conditions allow for the development of particular fauna and flora (Amorim et al. 2019; Lewis et al. 2019; Woolway et al. 2019; Amorim 2020; Bassone-Quashie et al. 2023).

Because of their social purpose, they are often found near urban and metropolitan areas, which expose them to the anthropic pressures typical of those areas. The scenario of climate change altering the atmospheric driving forces on a lake's hydrodynamics presents a new concern regarding how to predict impacts on the lake and understand what can be done to preserve its environment and social functions (Lewis et al. 2019; Woolway et al. 2019; Amorim 2020; Piccioni et al. 2021; Plec et al. 2021; Bassone-Quashie et al. 2023).

In this context, inland waters play a key role in climate change studies. This is especially the case with shallow water bodies because of their faster response to external force variability. Thus, the research on those environments helps enhance knowledge about climate change in all aquatic ecosystems (Amorim et al. 2019; Woolway et al. 2019; Piccioni et al. 2021; Yang et al. 2022).

Over time, researchers have sought to comprehend the behaviour and relationships of lake ecosystems using relational curves. These include air and water surface temperature (Toffolon et al. 2014; Piccolroaz et al. 2018), water temperature and aquatic life (Esteves 1988), specific mass, salinity, and water temperature (Boehrer & Schultze 2008), wind velocity and shear stress (Smith 1979), chlorophyll-a and phosphorus (Horne & Goldman 1994), chlorophyll-a and Secchi depth (Carlson 1977), phosphorus, detention time, and Secchi depth (Rast & Lee 1978), and a lake's stability and depth (Ambrosetti & Barbanti 2002). While these relational curves are effective, there are instances in which multiple relationships must be assessed simultaneously. In such cases, mathematical models became more prominent.

Mathematical models are nowadays a popular tool for studying climate change impacts in lakes. They are powerful tools that can be used to foster the better understanding, planning, and creation of forecasts of environmental systems (Bruce et al. 2018; Polli & Bleninger 2019; Amorim 2020; Piccioni et al. 2021; Plec et al. 2021; McBean et al. 2022; Papadimos et al. 2022; Sisay Kebede Balcha et al. 2023).

To obtain dependable outcomes, analysis tools, such as models or curves, require access to local data, comprehension of their underlying processes, familiarity with the tool, time to perform calibration and validation procedures, and the necessary knowledge to select the appropriate tool for each specific situation (Bruce et al. 2018; Polli & Bleninger 2019; Amorim 2020; Piccioni et al. 2021; Plec et al. 2021).

Choosing the appropriate tool to predict how climate change impacts water resources can be a challenge due to the various input data, set-up requirements, calibration coefficients, and types of results that each tool provides. To address this matter, this paper aims to compare the applicability of different tools used to predict climate change impacts on a lake's thermal stability.

Three tools were studied: an expeditious one, a medium one, and a slow one. The agility of a tool is classified based on the time it takes to set up and calibrate. The expeditious option is a relational curve that does not require calibration and can be set up with minimal input data. The medium option is a unidimensional mathematical model that automatically calibrates and requires only a few input data to set up. The slow option is a quasi-3D mathematical model that needs a large amount of input data and takes a long time to set up and calibrate.

Three distinct tools were used to analyse and evaluate the impact of two climate change scenarios on the hydrodynamics of the Hedberg reservoir: the thermal stability curve (TSC), a one-dimensional hydrodynamic model (General Lake Model), and a quasi-3D hydrodynamic model (Delft3D).

The three tools were applied over the studied area, using the same input dataset. The simulations’ results were compared, and the tools' applications were discussed in terms of quantity of input data; execution time; computational requisites; results' accuracy; and environment representation.

The lake has a 0.26 km² surface area, an approximate volume of 1.5 hm³, and a detention period of 2–10 days. Its maximum depth is 5 m and its mean depth is 4.5 m. In relation to its thermal regime, the Hedberg reservoir presents polymictic behaviour, with several stratification and mixing events observed throughout the year (Amorim 2020).

The lake's monitoring system is composed of two buoys, attached to a rope and a plummet, each one holding a vertical set of thermistors that gather temperature data every 2 min. One of the buoys is located at the entrance of the reservoir, near the main inflow (the Ipanema River), and the other is placed at the centre of the lake, where the depth is 4.5 m.

A local meteorological station is located on the lake's northeast bank and records air temperature, water level, atmospheric pressure, relative humidity, and precipitation every 10 min and solar radiation, wind velocity, and wind direction every 5 min.

To perform the climate change analysis, two scenarios proposed by the Fifth Assessment Report (AR5), published by the Intergovernmental Panel on Climate Change (IPCC), were applied. They are RCP 4.5 and RCP 8.5.

RCP 4.5, which predicts CO₂ equivalent emissions of 650 ppm and a mean temperature rise of 1.8 °C by 2100, was determined as the optimistic scenario. RCP 8.5 was delimited as the pessimistic one, predicting CO₂ equivalent emissions of 1,000 ppm and a mean temperature rise of 3.6 °C for the same period (IPCC 2014; Hölbig et al. 2018).

The climate change data were obtained from the PROJETA platform (Hölbig et al. 2018), a Brazilian government research platform responsible for assessing and providing access to climate change projections for South America. For South America, Central America, and Caribe, the scenarios RCP 4.5 and RCP 8.5 were generated with the regional climate model Eta, with a 20 km spatial resolution. The downscaling was forced by two global climate models, HadGEM2-ES and MIROC5 (Hölbig et al. 2018).

The first tool applied was the TSC. It is a correlation between atmospheric forces and water temperature that can determine the thermal status of the lake and its stability condition. The correlation was built based on the water column energy balance, analysing the energy retained in the water column (proxy by the Schmidt Number) with the incoming energy provided by the radiation (stratification energy) and with the outcoming by the wind's action (mixing energy) (Imberger 1998; MacIntyre & Melack 2009; Amorim et al. 2019; Amorim 2020).

One of the objectives of this tool is that it be based on a few easy-to-measure variables. The input data are incident radiation, water temperature, and wind speed. Incident radiation was used to calculate the sum of the radiation multiplied by the photoperiod, in the last 24 h time window, and it gives the specific stratification energy (J m⁻²) provided (Amorim 2020).

To calculate the mixing energy, the fluctuation of wind speed was used. Based on the concept that a constant velocity only affects the surface of the lake, while the movement induced by variations changes the way energy is transferred to the water, causing the water column to tilt. Thus, the mixing specific energy (J m⁻²) was proposed as the energy associated with the fluctuation of the wind velocity (Kullemberg 1976; Imberger 1998; Amorim et al. 2019; Amorim 2020).

After this assumption, the energy associated with the instantaneous part was determined by the mean of the wind's speed variance in a time window (last 24 h) (m s⁻¹), multiplied by the air density (kg m⁻³), and the lake's depth (m) (Amorim 2020).

The role of each intervening variable was represented by a proxy, namely:

• Rad: The total amount of the net incoming radiation at the surface in the last 24 h (J m⁻²) (Equation (1)).

Lakes' thermal condition is a result from which is the ruling energy of the balance (stratification or mixing). The proxies reflect the effectiveness of wind or radiation energy transfer to the water column, and so, which is the ruling energy on balance at the moment. They were correlated using the dimensionless ratios S* Rad⁻¹ and W* S*⁻¹ (Amorim 2020).

The one-dimensional hydrodynamic model selected was the GLM (version 3.0.5), developed by the Aquatic EcoDynamic group of the University of Western Australia to support the Global Lake Ecological Observatory Network (GLEON) initiative (Hipsey et al. 2019).

Extensive literature can be found on the GLM, as it is an open-source model widely used in scientific studies (Bruce et al. 2018; Hipsey et al. 2019; Woolway 2019). In the past, many works have applied the GLM to perform hydrodynamic simulation of lakes and reservoirs facing the impact of climate change scenarios (Huang et al. 2017; Farrell et al. 2020; Ladwig et al. 2021).

The GLM uses a deterministic, mechanistic, time-dependent, and numerical solving approach, through a Lagrangian structure, to simulate hydrodynamic processes that take place within the vertical dimension of the lake. The model can simulate the occurrence of stratification and mixing events. For that, the system's energy and mass balances are solved at daily and/or hourly timeframes, depending on the available input data. At every timestep, the vertical density profile is evaluated, as the Kinect energy (KE) inputted into the system is compared to the KE required to break the stratified structure.

The GLM allows for the calibration of various parameters to enhance the accuracy of the simulation. User-defined settings include maximum and minimum layer thicknesses, as well as maximum layer volume. Sources and sinks of water can be calibrated with scaling factors, inflow and outflow characteristics, the seepage coefficient, and the runoff threshold. For the energy balance, parameters such as the light extinction coefficient, sensible and latent heat transfer, the wind drag coefficient, and mixing processes efficiency coefficients are adjusted for the surface layer. Additionally, calibration parameters for the bottom layer encompass the diffusivity coefficient and the soil-sediment thermal conductivity.

In this study, the simulation timeframe used in the 1D model was from January 2017 to December 2020, with hourly time steps. The lake's initial depth was set at 4.1 m, as observed in the monitoring data, and the system salinity was set at 0.0005 mg L⁻¹. The GLM optional modes were rad_mode = 1 and deep_mixing = 1. Other set-up parameters were minimum layer volume (0.01 m³); maximum number of layers (900); and minimum and maximum layer thicknesses (0.01 and 0.7 m). Table 1 summarizes the boundary conditions applied in the model and the parameters to be calibrated. Parameters not specifically mentioned were set to the model default values.

The third tested tool was the quasi-3D hydrodynamic model performed in the Delft3D. It is a mathematical model based on the resolution of the Navier–Stokes equations, using the finite difference method. The main reasons for this choice were the model capacity of quasi-tri-dimensional modelling, which applies to the research objective, the widespread and known reliability of the software, and the fact that it is an open-source system (DELTARES 2014; Polli & Bleninger 2019; Plec et al. 2021).

The study site was described in the model as an orthogonal grid with 13 m cells, representing the spatial variations in the surface area, and the water column was described in 30 layers, with 0.2 m of thickness, in the z-model.

The boundary conditions were defined using the collected data. The physical and hydraulic parameters used were the lake's bathymetry, the upstream inflow, and the spillway rating curve. In the heat models, variables such as radiation, wind direction and velocity, air temperature, precipitation, humidity, and evaporation were used as the models' driving forces (Figure 2). The available variables and previous experiences (Soulignac et al. 2017; Polli & Bleninger 2019; Plec et al. 2021) justified the choice of Ocean as the heat model on the Delft3D.

Parameters must be set up in the model's calibration to adjust the simulated process to local characteristics. As presented in the 1D model, Table 2 summarizes the boundary conditions and calibration parameters used in the 3D model. Parameters not mentioned were set to the default values.

To perform the climate change simulation, the data from PROJETA were applied to the model boundary conditions, along with the rating curve of the spillway and the morphometry data of the lake.

For the above equations, n points to the sample size, e is the error, is the predicted value of the criterion, is the measured value of the criterion, and is the mean of the measured values.

The RMSE has been used as a standard statistical metric to measure model performance in meteorology, air quality, and climate research studies. Another useful and widely used coefficient in model evaluations is the MAE. The difference between them is that the MAE gives the same weight to all errors, while the RMSE penalizes variance, giving errors with larger absolute values more weight than errors with smaller absolute values (Chai & Draxler 2014).

NMAE is normalized to the mean, enabling like comparisons between variables, and is absolute, so that under and over-estimations do not cancel each other out (Bruce et al. 2018). Recognizing the limitations of the correlation coefficients, Nash & Sutcliffe, in (1971), proposed an alternative goodness-of-fit index, which is often referred to as the efficiency index (E_f). The advantage of the Nash–Sutcliffe index is that it can be applied to a variety of model types. For linear models, the efficiency index will be in the interval from 0 to +1. For biased models, the efficiency index may actually be algebraically negative. In nonlinear models, negative efficiencies can result even when the model is unbiased (McCuen et al. 2006).

In addition to the RMSE and MAE, Nash–Sutcliffe can be a useful index. However, it can be sensitive to a number of factors, including sample size, outliers, magnitude bias, and time-offset bias. Therefore, the literature suggests using a combination of the indices to assess the model, as was done in the research (McCuen et al. 2006; Chai & Draxler 2014; Bruce et al. 2018).

The hydrodynamic models were calibrated and validated so that they could be used in the climate change simulations. The periods used in this step are listed in Table 3. As mentioned earlier, gaps in climate data were filled by data from a nearby meteorological station and from the water level with a hydrological model developed for this watershed (Magalhães 2023).

The models were calibrated by applying the values presented in Table 4, and the results were both graphically and statistically evaluated. The models' accuracy indices MAE, Nash–Sutcliffe, RMSE, and NMAE (Table 5) showed good agreement between field measures and simulated results when compared to literature values (Soulignac et al. 2017; Polli & Bleninger 2019; Plec et al. 2021). For example, a study that modelled 34 lakes has an NMAE mean of 0.11, a maximum of 0.25, and a minimum of 0.04 (Bruce et al. 2018).

The results of climate change simulations for the three tools were consistent, indicating a reduction in the number of overturning events, an intensification of stratification, an increase in the number of days in this condition, and an increase in the energy required to break the density gradient created, thereby altering the thermal patterns once observed. Those impacts were also found in other studies regarding climate changes in lentic environments (Piccioni et al. 2021; Yang et al. 2022; Bassone-Quashie et al. 2023).

The gain that comes with the model's utilization is the capacity to represent the water column behaviour through the years and confirm the trend observed in the curve, i.e., that climate change will increase stratification events and the mean temperature of the lake. Also, with the three-dimensional model, it is possible to observe the surface spatial variability of the temperature, indicating the points most affected by the higher temperatures' impacts.

Comparing the results between the two models, it is noted that the three-dimensional model shows a warmer layer bottom than the unidimensional model in the pessimistic scenario. This is an important difference because, in a longer-term evaluation, the heating of the bottom layer can lead to a well-mixed water column throughout the year, changing the actual lake's mixing regime, which is polymictic.

One way to symbolize this change in the mixing regime is by quantification of the number of mixing events that occurred along the simulated period. This proxy was verified for the first 5 years of the climate change time series, considering that the lake's water column mixed when the difference between the bottom and the surface was less than 1 °C. The results showed fewer mixing events than the current situation, confirming the change in the lake's mixing regime an effect of the water column temperature's elevation.

A mixing regime with longer and more energy-accumulated stratification periods damages the gas exchange between layers and enhances dissolved oxygen depletion in the deeper layers. This condition makes the mixing events more harmful to the water quality. Each time the water column overturns, the deeper layers, with a higher load of organic matter and lower levels of dissolved oxygen, will be brought to the surface, impairing the water quality throughout the whole water column.

The three tools produced consistent and similar results, although their applicability might vary based on analysis objective, available time, and data. The purpose of this paper is to discuss this applicability, which is done in the next paragraphs.

The stability curve is a simple tool that requires minimum input data, uncomplicated equipment, and little computing time. Its results can help managers predict the trend of environmental stability many years ahead (up to more than one hundred years). This tool is considered lumped because it does not require any calibration and has a simple set-up process, which makes it less time-consuming than more complex tools.

The unidimensional model is simpler and requires less input data, time and specific knowledge to set up compared to 3D models. It can be calibrated automatically and has a short simulation time (only 1 min to simulate 5 years). Additionally, it provides an accurate representation of the water column, especially the water surface.

The three-dimensional model offers significant improvements in representing the water column and spatial variability. However, it requires a larger amount of input data and a more advanced understanding of the software. Setting it up can be time-consuming, as it involves building a grid with multiple boundaries, defining initial conditions, and determining where to print the results. Additionally, calibration is more complex because it involves verifying the behaviour of variables both horizontally and vertically and adjusting several parameters. Although simulations can take longer (e.g., 7 days to simulate 5 years), there are ways to speed up the process.

This analysis highlights the fact that forecast tools, when appropriately built and calibrated, can accurately represent lakes' hydrodynamics and predict the climate change impacts on them. Also, this analysis brings to light the differences between the applications of simpler and complex instruments.

In terms of choosing the correct tool to perform forecast scenarios, the comparison shows that this choice will depend mostly on the quantity and spatial variation of the data, the study objective, and the window simulation. A fast analysis with less detailed results can be provided with simpler data and tools, while a difference representation between upper and lower water column layers requires a minimum of a 1D model, which can be accurate only for small lakes with no major differences along their surface areas.

In the context of climate change and water scarcity worldwide, the development of tools to better understand, maintain, and improve water quality in lakes and reservoirs becomes an essential ally of environmental research and limnology.

The applications performed in this paper tested three different tools to forecast the influence of climate change on a lake's mixing regime and evaluated them in terms of input data quantity, the lake's hydrodynamic representation, and time consumption.

All methods presented an accurate and converging performance in representing the impacts of climate change on lakes' mixing regimes, showing the strong influence of the driving forces on the mixing regime and the warming of lakes. This situation leads to a reduction in the period of the colder season and promotes stratification.

Such effects will directly influence lakes' water quality, reducing the exchanges between layers, promoting anoxia problems on the bottom, and, consequently, making the mixing event problematic to the environment.

The pessimistic scenarios have mixing events with greater amplitude, which results from a longer and more energy-accumulated stratification in previous periods. This scenario creates higher vertical velocities, resuspending more organic load and dropping dissolved oxygen levels along the water column.

The tools' performance evaluations show that each has its pros and cons, which determine its applications. The TSC has the advantage of a faster response, a minor need for data input, and the ability to explore a larger simulation window. On the other hand, the models are capable of representing impacts along the water column.

The vertical unidimensional model requires less data and time than the three-dimensional model. Its limitation is the fact that it may accurately represent only small lakes with no major differences along their surface areas. If the study area has many changes in property along the surface and across the water column, a three-dimensional model is required, demanding more input data and time to perform the simulations.

Possibly, in lake management, it would be more appropriate to combine the methods, using the curve to analyse the trend faster and delimitate the exact period that requires more detailed study.

The authors would like to thank CAPES, USP, and FCTH for their support, data availability, and funding of research activities.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.