Analysis of different hypotheses for modeling air–water exchange and temperature evolution in a tropical reservoir

The construction, operation, and removal of large reservoirs are fundamental environmental issues because of the benefits they generate and the environmental impact they have. Hydroelectricity production accounts for 16% of total energy production (Killingtveit 2019) but, in the American continent, it is the largest source of energy production. In addition, they have essential functions such as supplying water to cities, irrigation, and energy. While these dams have brought some benefits, they have also brought severe and irreversible changes to the natural hydrology of the river, affecting soil, vegetation, biodiversity, landscape, atmosphere, and, in turn, watersheds (Agostinho et al. 2008; Kuriqi et al. 2021), reducing water flow, increasing water residence time, thermal stratification, increasing sedimentation rates and decreasing dissolved oxygen concentrations, among others (Pimenta et al. 2012).

The impact on reservoir water quality is a significant concern in ecological management. Life quality is highly constrained by water quality, which is affected by the quality of reservoirs built for hydroelectric projects and water supply. Building dams is particularly important for many Latin-American countries, where most of the electric supply comes from hydropower and many reservoirs. Thermal stratification caused by solar heating plays a significant role in determining water quality in the reservoir. In tropical countries, lakes and impoundments will stratify during the year as a result of increasing temperature differences between the warm upper (epilimnion) and cold lower (hypolimnion) layers of water (Halini Baharim et al. 2011). Air temperatures can vary by up to 10 °C from day to day, with water temperatures ranging from 20 to 32 °C.

Temperature studies in reservoirs have ranged from field measurements to computer tools, such as modeling. In the case of tropical reservoirs, from a data analysis and statistical point of view, research focuses on (1) describing and characterizing thermal distribution by field observations, such as the characterization of the stratification pattern and affectation of these thermal regimes in response to climate variables; (2) reservoir biochemistry and aquatic biology and the relations between them.

In Labaj et al. (2018), a data logger recorded temperature profiles hourly and measured every 2 m at a point at a maximum depth of 30 m in four lakes in Ecuador. They observed that thermal stratification was significantly positively correlated with air temperature and negatively correlated with wind speed across all lakes, and the most critical finding was that the high-resolution data showed that stratification did not breakdown overnight. In Lewis (1996), a theoretical comparison of the characteristics of several tropical and temperate lakes was carried out. The main finding of this study was that the comparison of tropical and temperate lakes has excellent potential to demonstrate lake function in general and is typified by non-seasonal substantial variations superimposed on seasonal cycles in most cases. Elçi (2008) studied the effects of thermal stratification and mixing on reservoir water quality at one point in the reservoir over 14 days using a water quality meter. Multivariate analysis was carried out on a data matrix of seven variables. The results showed that air temperature, lagged wind speed, and humidity influence variations in water quality parameters. In Rangel et al. (2009), temperatures are measured with a portable digital meter, in two main climatological seasons: cool-dry season and warm-rainy season, at a central reservoir point, at the surface, 1.5, 2.5, and 4 m. The main finding of this study was that the thermal pattern strongly influenced the vertical distribution of the phytoplankton community. Pajares et al. (2017) observed that temperature and oxygen stratification shaped the distribution of picoplankton. In Huszar et al. (2006), the data set includes 192 aquatic systems sampled on seasonal bases for at least 1 year using average values. They found significant differences in nutrient–chlorophyll ratios and thermal profiles between tropical and temperate climates. Amorim et al. (2020) monitor ten reservoirs in Brazil to study cyanobacterial blooms. They demonstrated that omnivorous crustaceans and total dissolved phosphorus mainly influenced cyanobacterial biomass. Solar radiation, air temperature, mixing zone, and salinity also significantly explain biomass behavior. These studies highlight the importance of thermal stratification as one of the main factors affecting variation within the water column of tropical lakes. However, in evaluating the thermal structure of the reservoir through data analysis, it is difficult to measure and verify the quality of the data, and it is difficult to propose an analysis and management scenario.

From a numerical point of view, 1D (one-dimensional) lake models have been built (Samal et al. 2009; Katsev et al. 2010) to assess ecosystem health, giving daily information at one point, and at a depth interval of 0.5 m for 3 years daily. They observed that the changes in the stratification regime in these waterbodies affect the water quality and the ecosystem's health, primarily based on temperature and dissolved oxygen parameters. In Crowe et al. (2008), one monitoring station is used to examine the chemical composition of the water and estimate transport time scales in reservoirs; the main finding was that seasonal temperature variability affects biogeochemical cycling in lakes.Rueda et al. (2006), using 15 temperature and light profiles taken in the reservoir throughout 1 year, show that the temporal variations of mean residence times occur not only at seasonal time scales but also at shorter scales. Several 2D (two-dimensional) models have been built to characterize thermal stratification and assess water quality. In Lindenschmidt et al. (2019), with monthly measurements, showed the influence of climate change on water quality reservoirs. Mesquita et al. (2020), Basso et al. (2021), Azadi et al. (2019) modeled and analyzed different scenarios that could cause eutrophication processes using monthly measurements. Ziaie et al. (2019) and Chuo et al. (2019) found significant differences in nutrient entrance and the relationship with algae blooms with monthly measurements. However, most works have analyzed the water quality and thermal dynamics without considering 3D (three-dimensional) processes.

1D and 2D models have the advantage that they are much faster than 3D ones when performing numerical calculations; they allow long-term simulations and are easier to calibrate and validate as they depend on fewer parameters. For example, the 1D model revealed daily air and water temperature relationships but failed when the advective term was significant. It also fails to quantify the potential effects of climate change on water temperature in a shallow reservoir (Gooseff et al. 2005). Note that this approach cannot predict stratification dynamics. A 2D model is essential for studying stratification dynamics subject to horizontal advection. Even though 1D and 2D models can consider short- and long-wave radiation in the thermal model, heat flux by evaporation/condensation, and the convection process at the free surface, they cannot capture the processes of heat exchanges with both the bottom and atmosphere, while the 3D does. Besides, the turbulence models have been demonstrated to be important in maintaining the thermoclines and are fully developed only for 3D models (Goudsmit et al. 2002).

In the case of 3D models, phenomena at multiple scales and couplings can be modeled. Other 3D flows, including gravity current, significantly affect the water quality and temperature (Kopmann & Markofsky 2000). In complex geometries with highly three-dimensional flows, the more traditional depth or width-averaged models cannot accurately capture mechanisms affecting temperature transport and mixing (Politano et al. 2008). Overall, 3D modeling tools are critical to understanding the interaction between all aquatic ecosystem components. 3D modeling enables a better assessment of the high complexity of the reservoir and its natural cycles.

Regarding 3D modeling, there are several studies for shallow reservoirs, such as study of the influence of cold fronts on the heat fluxes and thermal structure in a tropical reservoir (Curtarelli et al. 2013), study of the circulation patterns in a shallow tropical reservoir (Yang et al. 2019), study of saline intrusion (Laval et al. 2005). However, while numerous studies have been conducted to characterize thermal stratification in temperate lakes and oceans, there have been fewer numerical studies of temperature dynamics in tropical reservoirs (Politano et al. 2008).

3D modeling has proven to be a valuable tool to advance the understanding of the physical processes of fluid dynamics and water quality, thermal processes such as the effects of wind and temperature induced flow (Matta et al. 2017), the effect of the wind and tidal forcing (Moloney et al. 2016), thermal discharge released to coastal areas (Gaeta et al. 2015, 2020), coupling with complex water quality process (Piccioni et al. 2021). It has also been used to assess environmental impacts from cooling and heating power plant production (Ligier & Okumura 2019), to analyze flood and ebb events and analyze indicators of the contamination degree (da Silva et al. 2021), to evaluate Escherichia coli development scenarios (Bedri et al. 2013), to study the dynamics of geomorphological features (Lisboa & Fernandes 2015), or the heating impact of a reservoir on downstream water temperature (Jiang et al. 2018; Plec et al. 2021). However, all these simulations have been performed without heat exchange with the atmosphere and not for tropical reservoirs.

Other 3D models use a simplified heat exchange model (Scanlon et al. 2020) to assess the climate change scenario's impact over a small tropical lake (Duarte et al. 2021), to evaluate the impact in coastal areas, simulate flooding and hydrodynamic patterns in wetlands (Costi et al. 2019), and to study thermosaline circulation in an estuary (Vouk et al. 2019). However, as explained by Politano et al. (2008), these studies are not performed on tropical reservoirs, and other 3D models have been developed to use the entire atmosphere–water exchange considering hydrometeorological conditions (Angelotti et al. 2021) but without complex geometry.

Some recent studies have assessed the effects of the advection schemes and turbulence and other important numerical parameters with salinity as a tracer on 3D models (Smolders et al. 2015; Chen 2020; Justin-Brochet et al. 2021), to evaluate the non-hydrostatic 3D model in rivers (Politano et al. 2008), optimize numerical parameters in 3D models for the calibration in small lakes (Merkel 2019), or coastal areas (Cooper & Spearman 2017; Gaeta et al. 2020), but they were not used to analyze those conditions in a stratified deep tropical large reservoir with spatiotemporal data of high frequency for calibration.

Complex reservoirs such as deep tropical lakes, with elongated heterogeneous basins, with the influence of several rivers or advective flows (Marín-Ramírez et al. 2020), are vertically complex systems, making difficult both the building and the implementation of a model that correctly represents the thermal structure.

In this sense, the present work emerges as a response to integrating a robust water quality model that studies the effect of hydroclimatological variables such as air temperature, wind, and internal mixing process on the thermal dynamics of a tropical reservoir in Colombia, using highly resolved simulations of hydrodynamic, thermal, and energy processes. The importance of these processes has been demonstrated by data analysis from the observations of La Miel reservoir (Alzate-Gómez et al. 2023).

Based on the context discussed here, the study analyses thermal dynamics over the vertical water column in a tropical reservoir. The general objective of the present work paper is to show the potential and limitations of three-dimensional hydraulic models in predicting the thermal distribution of rivers with reservoirs and hydro-plants. This leads to analyze different modeling hypotheses for temperature and air–water exchange and evaluate the model performances on short periods according to the observations. Section 2 outlines the characteristics of the study site. Then, Section 3.1 presents the available data. Section 3.2 describes the governing equations, and Section 3.3 details the model that has been implemented along with the performance criteria. Finally, Section 4 describes the results, and a discussion of the main findings is given in Section 5.

‘La Miel I’ dam is a gravity dam located on the riverbed (see Figure 1) and was built between 1997 and 2002 for hydroelectric power generation. The tropical reservoir Amaní (5.25° lat, −75° lon) is upstream of the dam. The dam is 188 m high, giving it a storage capacity of 571 million m³ and a surface area of 12.2 km² with a generation installed capacity of 396 MW in three turbine units. Its commercial operation began in December 2002. In 2010, the Guarino diversion dam on the Guarino River was opened, and the Manso diversion dam on the Manso River began operations in 2013. Both divert water into the ‘la Miel I’ Reservoir through a tunnel.

Analyzing the data from the meteorological measuring stations is described in Section 3.1 ‘Field Data and the Monitoring System’: the average temperature is 25.3 °C, the maximum temperature is 33 °C, and the minimum is 18 °C. The highest monthly evaporation rates occur in the dry season, with 188 mm/day, while the lowest evaporation rates occur in the wet season, with an average of 56 mm/day. Relative humidity varied between 40 and 100%, with a mean value of 80%.

The annual cycle of rainfall in the area presents a bimodal behavior, with maximum flows in November–December and March–April and recession in July–August and January–February, with the recession in July–August being considerably more substantial. The lowest precipitation values occurred in July at 122 mm and the highest in October and November with values of 416 mm. The average flows of the La Miel and Moro rivers at the Puente Samaná and Tarro Pintado (see Figure 1 stream points) gauging stations are close to 50 and 20 m³/s, respectively (estimates made with data from 2003 to 2016).

In this work, high-frequency field data were used, 3D simulations were performed, and the former was compared with the latter. The data needed for the implementation of the model were analyzed, the mathematical equations and their implications for the modeling were analyzed, and the model configuration was adjusted accordingly as a basis for the understanding of the thermal dynamics. The hydrodynamic model was calibrated according to the observed data and the error was quantified using metrics. Then, the model responses to different thermal forcing are analyzed and the results were compared with the observed data.

The materials and methods used for this purpose and the specifications for each stage are explained in Sections 3.1, 3.2, and 3.3.

The monitoring system comprises water quality measurements at five stations, of which only the E7 station next to the dam will be used, in addition to thermistor chains and CTD (Conductivity Temperature Depth) instruments. For the same period, five hydrological parameters (flow rate, precipitation, evaporation, wind, and air temperature) covering all areas of the basin have been collected; specifically, seven monitoring stations have been used to measure precipitation, two stations for flow rate and two stations for the other hydrometeorological parameters (gauging stations, Figure 1).

In station E7 as shown Figure 2, water depth with ADCP (Acoustic Doppler Current Profiler) instrument and temperature measurements were made at five different depths: two in the epilimnion (subsurface and center), one in the middle zone, and two in the hypolimnion (center and three meters from the bottom). In addition, a thermistor (Figure 2) was installed near the dam (E7) to have a continuous temperature record from April to September 2016. In addition, four thermistors were placed in this reservoir at 1.5, 4.5, 15, and 40 m depth and allow to measure the typical behavior of the temperature profile. The thermistor located at 1.5 m records the temperature of the superficial mixed layer (epilimnion), the one at 4.5 m is located approximately in the thermocline, and the one at 15 m records the temperature of the hypolimnion above the intake structure in the layer that presents the most significant temperature variations and the 40 m to capture the temperature at the bottom. All thermistors were programed to collect data every 15 min. These measurements were taken between May 2016 and September 2016.

Additionally, three field campaigns were carried out to collect temperature on about 36 profiles (CTD points Figure 2) in the reservoir and tributaries. The campaigns were carried out on the following dates: (1) 6 and 7 July, 2015 (monthly monitoring), (2) 16 and 17 August, 2015, and (3) 22 September, 2015. Vertical measurement profiles were taken at different reservoir points using a master CTD RBR profiler with temperature, turbidity, conductivity, and chlorophyll sensors. This equipment allows data to be taken with a frequency of 6 Hz along the vertical profile line. The obtained profile shows a resolution of about 5 cm.

Hydrological parameters were collected hourly for 183 days from 20 May, 2016 to 30 September, 2016 (El Vergel gauging station, Figure 1). Hydrometeorological data were obtained from the DHIME website (http://dhime.ideam.gov.co/atencionciudadano/) of the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM, acronym in Spanish). The meteorological information used for the model was provided from the gauging station closest to the dam.

Some missing data as flow discharge from the other tributaries, Santa Bárbara, San Luis, and Fraguas, and water temperature for all the tributaries, were calculated from the SWAT (Soil and Water Assessment Tool (Gassman et al. 2007)) model (see Figure 2 stream points). SWAT is a semi-distributed and physically based model for daily simulation of daily discharge, sediments, nutrients, and water quality parameters. The most critical processes simulated by the model are surface runoff, infiltration, lateral flow, base flow, evapotranspiration, and groundwater recharge (Tessema et al. 2021). The model was calibrated using daily discharges from the Puente Samaná and Tarro Pintado gauging stations from 1985 to 1990 and validated from 1990 to 1997.

In Figure 4, the wind speed data are presented. During the study period, the minimum wind was 0.05 m/s, the average wind was 0.73 m/s, and the maximum was 3.0 m/s. The canyon topography of the area reservoir may explain the low wind speed values. The reservoir operated at a low water level in August, and the water level varied by 7 m, reaching 406.5 m on 27 August. The water level remained like this for 2 weeks, then rose again until it reached 417 m by the end of September (Figure 5). The inflow remained roughly constant during August, at 25 m³/s into Miel, except at the beginning of September, where the peak flow was as high as 65 m³/s in La Miel's arm. The reservoir presents indeed more inflow from the La Miel arm, as seen in Figure 6.

The water temperature distribution in the analysis process presents no mixing phenomenon, and there are at least six layers in August–September. The layers that present the most significant variation are in the first meters up to 10 m. The temperature distribution of the first four layers varies greatly, while the temperature distribution of the two deepest layers does not change much. Also, the high-resolution data suggest that stratification does not breakdown overnight.

The vertical structure of water temperature is stable in different periods of each month, and the surface temperature is relatively high throughout the year, with a mean value of 29.6 °C and a maximum of 31.8 °C at the free surface. At depths between 25 and 15 m, the temperature in the reservoir appears moderated by the temperature entering the tributaries. The temperature in the deep layers is between 22 and 23 °C. The analysis in the next section indicates that the external heat flux, which is determined by meteorological and inflow conditions, is insufficient for altering the nearly stratified state of the reservoir.

Clearly, there is a permanent stratification in the reservoir: hourly and daily fluctuations of the isotherms 28, 29, and 30 °C in the epilimnion (of about 1–3 m) are observed. On the monthly time scale, we also observe an overall increase of the temperature of about 2 °C in the dam followed by a decrease of similar intensity, as the isotherms plunge by about 5 m and get back to their initial location.

One of the critical factors in the formation of the thermal structure is the surface heat flux, which is composed of solar radiation, incident and emitted long-wave radiation, sensible heat flux, and latent heat flux. Often, this energy flow results in an increase in the temperature of the surface layer of the reservoir and, thus, a decrease in its density. In contrast, the deeper layers remain cooler because they are less exposed to these heat sources as shortwave radiation penetration is limited, and diffusive heat transport is too slow, as seen in the measurements (Figure 8). In the tropics, the absence of seasons means that the difference between maximum and minimum annual radiation is slight. In addition, the total daily radiation is remarkably constant throughout the year since there are no significant variations in the number of daylight hours per day. This lack of variation in the daily radiation received means that surface heat fluxes do not vary so dramatically, and therefore, monthly patterns of stratification and mixing are not very pronounced in the present case.

From August to September, some variations in the reservoir temperature appear due to changes in the inflow and outflow of water, especially the difference in the operation level of the reservoir and meteorological conditions. Hence, the water temperature stratification pattern slightly differs from August to September as shown Figure 8.

In order to achieve the objectives outlined, it is necessary to have a 3D model of coupled hydrodynamics and temperature, which must be validated to ensure that what is simulated is representative compared to observed behavior.

For the hydrodynamic component, the code solves three-dimensional hydrodynamics equations with a coupled energy equation regarding temperature. The equations presented were obtained with the following assumptions (Hervouet 2007):

• - Three-dimensional equations of incompressible flow with a free surface changing in time (Equations (1) and (2)),

• - Negligible variation of density in the conservation of mass equation (incompressible fluid),

• - Hydrostatic pressure (Equation (3)),

• - the Boussinesq approximation (Equation (6)) for the momentum (the density variations are only considered in the buoyant forces).

Here, U (m/s) indicates the three-dimensional component of velocity, t (s) indicates the time, (kg/m³) indicates the density, (kg/m³) indicates the reference density, g (m/s²) indicates the acceleration due to gravity, p (Pa) indicates the pressure, F (m/s²) represents the external forces, h (m) indicates the water depth, Z_S (m) indicates the free surface elevation, p_atm (Pa) indicates the atmospheric pressure, (m²/s) indicates the eddy viscosity, x, y (m) indicate the horizontal space components, z (m) indicates the vertical space component.

The parameter b depends on the location site and will be calibrated in this work, V_wind is the wind velocity (in m/s), and T is the water temperature at the free surface in °C.

The unstructured mesh comprises 4.779.814 nodes and 8.909.098 triangular elements with a minimum size of 0.7 m in a 2778.4 km² area, created with the software Blue Kenue. The mesh is made of prisms representing the spatial variations in the surface area. The model mesh is developed as a series of planes between the bed and the free surface. Twenty-three planes are used to describe the water column. These layers have a percentage distribution of 10 m at the bottom to 0.5 m near the free surface. It is a classical sigma transformation that follows the shape of the terrain (TELEMAC3D 2020).

As boundary conditions this research imposes a constant temperature at the inflows (of T₁ = 22.9 °C, T₂ = 22.6 °C, T₃ = 24 °C, T₄ = 23.4 °C, and T₅ = 23 °C) (see Figure 10). The boundary condition imposed on the free surface was Dirichlet for one case study and Newman for the other (Equation (12)), using air temperature, wind direction, and velocity based on data collected from the local and nearby weather stations. These approaches are developed in detail in Section 4.2.

The modules employed were TELEMAC-3D-WAQTEL version v8p2 (TELEMAC-3D 2020). The simulations were conducted for each case investigated through parallel computing on the supercomputer OLYMPE from the CALMIP computing center for public research in Toulouse. The numerical parameters were defined by trial-and-error, where the model was run several times and yielded the following results. The set of selected numerical and physical parameters is summarized in Table 1.

To solve all processes, including hydrodynamics and thermals, the GMRES solver was chosen based on mass conservation and computational cost. To simulate all the physical processes, this study uses a constant viscosity model for the turbulence model (Hinkelmann 2005; Hervouet 2007). The turbulent viscosity value was fixed at 10⁻³ m²/s in all directions, as a results of previous simulations using Smagorisky as a turbulence model. Several simulations were run for a time step ranging from 0.1 to 5 s. The time step was set to 1 s since it provided the best compromise between stability and performance.

The Multidimensional Upwind Residual Distribution (MURD) method was applied for the advection of three-dimensional variables under TELEMAC-3D, and the boundary conditions were applied following the method of characteristics. The advection of tracers was solved using the distributive MURD N method. Section 5 discusses the parameter choice.

In the equations, n is the sample size, is the predicted value of the output (water level or temperature) at time or location i, X^o_i is the corresponding measured value, and is the mean of the measured values.

The MAE and RMSE have been used as standard statistical metrics to measure model performance in meteorology, air quality, and climate research studies. The difference is that the MAE gives equal weight to all errors. At the same time, the RMSE penalizes the variance, giving more weight to errors with larger absolute values than to errors with smaller absolute values. As a result, when both metrics are calculated, the RMSE is never lower than the MAE (Chai & Draxler 2014).

Another critical aspect of the error metrics used for model evaluations is their capability to discriminate among model results. In this case, the advantage of the Nash–Sutcliffe index is that it can be applied to a variety of model types. This efficiency index ranges from −∞ to +1, 1 being the best. As well as RMSE and MAE, the Nash–Sutcliffe is a helpful index; however, it can be sensitive to several factors, including sample size, outliers, magnitude bias, and time-offset bias. So, using a combination of those indexes to assess the model is better. These statistics are commonly used to evaluate the model's goodness that fits the observed data (Lindim et al. 2011; Amorim et al. 2021).

Considering that the model is very sensitive to the temperature forcing given by the atmosphere–water exchange, it is necessary to analyze the thermodynamic behavior of the exchange process concerning the different forcings. Based on this, the following approaches are defined for the temperature boundary condition at the free surface:

• A.

  A constant water temperature T_water, of 29.8 °C, at the free surface (Dirichlet) without exchange with the atmosphere, considering the stability of the thermal stratification according to field measurements. The approach A was simulated until a steady state was reached. This approach served as initial conditions and reference for comparison with the other approaches.

• B.

  Exchange with the atmosphere at the free surface boundary condition (Equation (12)) with meteorological forcing (air temperature T_air and wind velocity V_wind) constant in time. A value of 32 °C is imposed for T_air based on the maximum value measured for the day of the analysis. V_wind = 0.6 m/s refers to the mean value over the period.

• C.

  Exchange with the atmosphere at the free surface boundary condition (Equation (12)) with meteorological forcing changing over time. This last approach uses observed hourly data of T_air and V_wind over the period.

The proposed approaches simulate the temperature profile in the reservoir under different forcings and parameters that best represent tropical conditions. They will be compared with the field thermistor data covering the period from 11 August to 19 August, 2016.

The same boundary conditions concerning the flow, source or intake, and tributaries temperature were set up for the three approaches; they have been detailed in Section 3.3.2 ‘Initial and boundary conditions’. These approaches were chosen for the sake of simplicity (number of parameters), data availability, and computation time.

The results of the temperature model for the three approaches mentioned in Section (3.3.5), which differ in the way they treat the water–atmosphere exchange, will be presented in Section 4.2.2. Considering the model's sensitivity to the parameters of this exchange, it is necessary to perform a detailed analysis of the impact of the physical parameters of the different models.

Sensitivity analysis of the thermic diffusivity of water (⁠⁠) and the coefficient of the atmosphere–water exchange model b is required according to the governing equations (Equations (12) and (13)) of the thermal model described in Section (3.2). Table 2 presents the values tested for these coefficients in each approach. For the local sensitivity analysis, the simulations were performed using the ‘one-at-the-time’ approach (Rajabi et al. 2015) by increasing each parameter by a given percentage while leaving all others constant and quantifying the change in model output.

Table 2

Simulations performed in this work

Case name .	Approach .	b .	(m2/s) .	Comment .
NO EXCH 1	(A)	N/A	1.0 × 10−6	The Prandtl number Pr = νt/κT = 1
NO EXCH 2	(A)	N/A	1.4 × 10−7	Pr = 7
NO EXCH 3	(A)	N/A	0	Pr = ∞
CAL0	(B)	0.0025	1.0 × 10−6
CAL1	(B)	0.0035	1.0 × 10−6
CAL2	(B)	0.01	1.0 × 10−6
CAL3	(B)	0.1	1.0 × 10−6
CAL14	(B)	0.01	0	Vwind = 0 m/s
CAL 1C	(C)	0.0035	1.0 × 10−6
CAL 2C	(C)	0.01	1.0 × 10−6
CAL 3C	(C)	0.1	1.0 × 10−6
CAL4	(C)	0.03	1.0 × 10−6
CAL5	(C)	0.05	1.0 × 10−6
CAL15	(C)	0.1	0	Taking as IC CAL3 last time step
CAL17	(C)	0.0025	1.0 × 10−6

Case name .	Approach .	b .	(m2/s) .	Comment .
NO EXCH 1	(A)	N/A	1.0 × 10−6	The Prandtl number Pr = νt/κT = 1
NO EXCH 2	(A)	N/A	1.4 × 10−7	Pr = 7
NO EXCH 3	(A)	N/A	0	Pr = ∞
CAL0	(B)	0.0025	1.0 × 10−6
CAL1	(B)	0.0035	1.0 × 10−6
CAL2	(B)	0.01	1.0 × 10−6
CAL3	(B)	0.1	1.0 × 10−6
CAL14	(B)	0.01	0	Vwind = 0 m/s
CAL 1C	(C)	0.0035	1.0 × 10−6
CAL 2C	(C)	0.01	1.0 × 10−6
CAL 3C	(C)	0.1	1.0 × 10−6
CAL4	(C)	0.03	1.0 × 10−6
CAL5	(C)	0.05	1.0 × 10−6
CAL15	(C)	0.1	0	Taking as IC CAL3 last time step
CAL17	(C)	0.0025	1.0 × 10−6

A, B, and C refer to the scenario type (see Section 3.3 for details).

View Large

A sensitivity analysis was performed on thermic drivers for air–water exchange included in the simulations to define which of them primarily influenced the models' performance as being sensitive to the variation of the simulated processes.

The analysis was carried out at the four depths of the dam, where Depth 1 corresponds to the temperature at the free surface, Depth 2 corresponds to the temperature at 4 m, Depth 3 corresponds to the temperature at 15 m, and Depth 4 corresponds to 40 m. These points were selected to represent the stratification; the field measurement data used to evaluate the simulations were taken between 11 August and 19 August, 2016. A spatial assessment compared to the temperature profiles measured by the CTD for 17 August was also performed. To investigate model performance, RE, RMSE, and MAE between measured and simulated temperature values were calculated.

According to Figure 17, the model could represent the main observed patterns, the first layer being the constant water temperature boundary condition at the free surface: (i) the second layer shows a slight observed variation in temperature between 28.0 and 31.41 °C for the period at an hourly scale. At 4 m after a transition period, NO EXCH 1 and NO EXCH 2 exhibit almost the same results with an increase of temperature which is not observed, whereas NO EXCH 3 simulates a constant temperature; the same behavior appears in the simulated results at 15 m; (ii) the cooling at the reservoir bottom, shows an overestimate of 0.5 °C for NO EXCH 3 and more than 2 °C and increasing with time for the other cases. For the cases of NO EXCH 1 and NO EXCH 2, there is a tendency for the bottom temperature to increase, which is not representative of the current situation for reservoirs with little change over time; (iii) the duration of the stratification period that never breaks in all cases; (iv) using approach A, it was not possible to reproduce the temperature oscillations, only the mean values which is due to the constant water temperature boundary condition at the free surface.

Unlike many natural lakes that experience overturning depending on weather conditions, no overturning was observed during the study at Amani Reservoir because it is a tropical lake, and the surface water temperature never drops below 16 °C. During the study period, strong stratification was observed in the Amaní reservoir, with a temperature difference between the surface and the bottom reaching 8 °C according to both measurements and simulations. The water temperature difference between the surface and bottom was 4 °C in NO EXCH 1 and NO EXCH 2 along the reservoir in all simulated times, while its maximum of 7 °C was reached at NO EXCH 3. The thermal structure of simulation NO EXCH 3 exhibit fewer layers of stratification than actually observed.

Approach B explores what happens to the temperature profile for different values of the exchange coefficient A (Equation (13)) with constant meteorological forcing. The heat transfer coefficients or A value reported in (Williams 1963) for reservoirs present a range between A∈[ 113, 378] W/m²/°C. It was not possible to find more recent values reported in the studies. As shown in Table 2, this study varied b (involved in Equation (13)) in the range [0.001, 0.1]. The latter values correspond to a range of A∈ [15, 900] W/m²/°C.

Generally, the model has the most difficulty reproducing the behavior of the second layer and the bottom (Depths 2 and 4). The best results according to the metrics are discussed in detail in the next paragraph.

Temporally, water temperature changes depend on the A air–water exchange coefficient. When A is lower, there is less exchange, and the water temperature experiences a drop, while when the A coefficient is higher, the water temperature rises. The case CAL2 is the one that is closest to the values observed at the free surface with an average of 28 °C; for the deeper layers, the A value that better represents the mean values is CAL 3.

Vertically, strong stratification was observed in the reservoir for those scenarios, with the surface-bottom temperature difference reaching 8 °C in all the cases according to Figure 22. In agreement with other studies, stronger stratification has been reported in other deep reservoirs. Among them, the temperature difference between the surface and the bottom of the reservoir reaches 10 °C (70 m deep), and the maximum temperature difference of the other reservoirs even reaches 14 °C (200 m deep) (Plec et al. 2021). One main reason for these gradients is the depth of the Amaní reservoir with 120 m. By preventing the penetration of solar radiation and wind-driven eddies, large bodies of water like the Amaní reservoir impede heat transfer and make it difficult to achieve a homogeneous temperature. These variables, included in the A coefficient, represent radiation, air convection in contact with water, and latent heat. Its correct determination improves the simulation results compared with approach A from a temporal and vertical point of view.

Longitudinally, the simulations of the temperature at the free surface drops at the entrance of the tributaries, and as the tributaries flowed toward the reservoir all the way down the dam, the simulated water surface temperature increases by 3.5 °C. Unfortunately, observed data do not provide any information on this phenomenon. The latter means that the heat transport in the reservoir is mainly due to advective transport, as was hypothesized, and in the first layer, it is due to the exchange of air–water and the atmospheric variables, mainly the exchange between air and water. Water is a poor conductor of heat, so thermal diffusion is slow in a reservoir. Therefore, the flow transfers energy with very high specific heat.

The cases analyzed for approach C reproduces the observed temperature oscillations at the free surface and the 4 m layer, but they differ in the standard deviation temperatures. Data measured at the free surface exhibited temperature oscillations of up to 1.5 °C, compared to 2 °C for the studied CAL1C case, 3.5 °C for CAL2C, and 4.5 °C for CAL3C. Indeed, the amplitude of daily temperature variation decreases with A. Temporally, water temperature changes with air temperature. When the air temperature was lower than the water temperature, the water temperature experienced a drop, while when the air temperature was higher than the water temperature, the water temperature rose.

According to Figure 25, the model could represent the main observed patterns: (i) the value at hourly time scale on the first layer of the surface, with a daily heat loss of about 4 °C for CAL1C, about 5 °C for CAL2C, and about 8 °C for CAL3C, for the period of analysis; (ii) at 4 m only CAL3C could reproduce water temperature oscillations but the average temperature is underestimated; (iii) cooling at the bottom of the reservoir revealed an overestimation of 0.4 °C for CAL1C and CAL2C and 0.7 °C for CAL3C.

In the case of approach C, the temperature oscillations at the free surface and at 4 m are simulated but the model still suffers from excessive heat losses, as the temperature of the layers below the air surface shows too significant daily decreases, which is not the case according to the observed data, as the air temperature shows oscillations of up to 14 °C. In contrast, the observed free surface temperature in the reservoir shows much smaller oscillations.

This study aims to understand the spatiotemporal exchange of tropical thermodynamics with the atmosphere. There are simulation tools that allow us to dig deeper into issues and knowledge, and potentially address scenarios such as climate change. However, when these tools are applied to an area as large as this case study, it often becomes a computational problem rather than a physical problem.

First, the use of a supercomputer includes decisions such as the number of processors and the desired computational time. This research used the OLYMPE supercomputer of CALMIP, supercomputing center which is providing resources and services for the entire scientific community of the Université Fédérale de Toulouse Midi-Pyrénées (UFTMiP). OLYMPE has a SEQUANA Cluster (ATOS-BULL) 1,365 Pflop/s, 374 computational nodes (36 cores/node), Intel^® Skylake 6140 Processor at 2.3 Ghz 18-cores.

For hydrodynamics simulations, 36 cores showed real-time performance with a ratio of actual computational time to simulated time of 0.8–1. However, the thermal equations solution required more cores to achieve similar performance when solved simultaneously with hydrodynamics.

Different configurations were tested using 36, 64, 120, 180, 288, and 324 cores. Using 324 processors makes simulations very fast, e.g., it took 16–18 h to run 8 days of thermic simulations; however, increasing numbers of cores exhibit diminishing returns, and this study believes that using 180 cores (which takes approximately 24 h to run an 8-day simulation) provides the best balance between runtime and computing resources.

From the hydrodynamic perspective, several factors were analyzed, and it was found that specific parameters did not affect the thermal distribution, such as the use of hydrostatic pressure, horizontal turbulence models (Constant Viscosity, Smagorinsky, k–ω, and k–ε), coefficient for horizontal diffusion of velocities and the friction coefficient. While the number of vertical planes, Boussinesq approximation, the boundary conditions on the free surface, the vertical turbulence models (constant viscosity, mixing length, Smagorinsky, k–ω, and k–ε), and the value of the coefficient for vertical diffusion of velocities influenced the temperature distribution in the vertical direction. In several cases, the mixing-length turbulence model results in faster mixing in the vertical direction than the observed thermal profile data, while constant viscosity best represents the thermal distribution in the vertical direction.

Regarding the thermal dynamics, this research found that the structure of the solution was very sensitive to the choice of an advection scheme. Seven methods for linear advection treatment (1) Method of characteristics, (2) Streamline Upwind Petrov Galerkin (semi-implicit), (3) Explicit finite volumes, (4) Explicit scheme + MURD (MURD) N scheme, (5) Explicit scheme + MURD PSI (TELEMAC 2020) (Positive Streamwise Invariant) scheme, (6) Explicit Leo Postma scheme for tidal flats, (7) Explicit scheme + MURD N scheme were used to analyze the behavior of the solution. Only the most relevant simulations were shown here for the sake of conciseness. Compared to other schemes, it was found that the MURD N scheme produces fewer numerical oscillations, less diffusion perpendicular to flow and thermodynamics, and is less computationally expensive. The scheme is also characterized by being conditional on a Courant number less than 1, satisfying the stability criterion.

The exchange process with the atmosphere is an essential and well-studied physical process. The tools used for its mathematical representation are still very limited, e.g., in the first layer of water temperature, the exchange between air and water also explains the behavior. The observations indicate that the external heat flux, which is determined by meteorological and inflow conditions, is insufficient to modify the nearly stratified state of the tropical reservoir. Many formulae in the studies can be used to compute this term (Equation (10)) but selecting an appropriate one is debatable. Generally, the net heat flux is calculated based on the water temperature and various meteorological data, primarily based on air temperature, air humidity, air pressure, measured shortwave solar radiation, and wind speed (although many others may be taken into account). An estimate of the net heat flux at the water–air interface turns out to be a challenge, especially when it comes to 3D modeling.

The first problem concerns obtaining the necessary input data for the analyzed site. Intuitively, the nearest meteorological station is the best source for meteorological data. Unfortunately, in many cases, the nearest station does not provide all the necessary data or is still too far from the considered study zone. Another problem is the differences in data obtained from the neighboring stations. Even in conducive situations when the necessary meteorological data are available near the reservoir, there are still uncertainties related to the station's location, for which conditions such as shading or wind speed are often considerably different from those at the reservoir.

The problem is widely discussed in the work (Benyahya et al. 2010). Moreover, some measured quantities may vary along the reservoir (or even across the reservoir width). Finally, very different values of net heat flux may be obtained. Williams (1963), for example, measured the heat fluxes for several reservoirs in different conditions. The results showed that the final sum of heat fluxes measured ranged between 113 and 378 W/m²/°C. In this case shows that A coefficient exchange value depends on the meteorological conditions, if they are constant over time or changing, and this research found the approximate values that can be used in tropical environment reservoirs.

Additionally, each term in Equation (12) is sensitive to several factors, and the chosen computation method, i.e., different formulae, may lead to varying results since they often depend on not well-defined parameters, factors, or coefficients site-specific. Moreover, in the ‘competition’ for the best formula, increasingly ‘more accurate’ formulae take into account more and more factors and thus require more and more input data, which, again in practical applications, are rarely available or highly uncertain.

In this work, different modeling assumptions for tropical reservoir temperature and air–water exchange are analyzed based on observations. This study evaluated the performance of various hypothesis in resolving hydrodynamic, thermal, and energetic processes, including the effect of hydroclimatological forcings such as air temperature, wind, and the internal mixing process (turbulence), and different boundary conditions in the free surface. The following conclusions can be drawn:

Reservoir water temperature is affected by heat exchange between the water and the atmosphere. This exchange is evident in the upper layer (O (<10 m)), and it is this dynamic in the first layer governs the distribution within the water body in areas where there is no advective flow from tributaries, which for the reservoir analyzed is on the order of 80%.

It was observed that the reservoir is permanently stratified and, unlike other reservoirs in the temperate zone, is not significantly affected by meteorological conditions outside the water body. This means that the thermal dynamics inside the reservoir remain very constant. In this work, numerical methods are used as a tool to study exchange. This research finds that a good definition of the boundary conditions at the water-free surface is crucial for a correct representation of the thermodynamics inside the water body.

The hydrodynamic TELEMAC-3D model has been satisfactorily calibrated using one month of daily water level observations at one point in the reservoir. Sensitivity analyses highlighted the importance of vertical turbulence models and vertical diffusion coefficient for a correct representation of the observed behavior.

Overall, reservoir water temperature is strongly influenced by heat exchange between water and the atmosphere. Accurate modeling of such phenomena is thus necessary. Three different approaches have been tested to study the impact of air–water exchanges at the free surface: a constant water temperature without exchange with the atmosphere (approach A), meteorological forcing with atmospheric parameters constants in time (approach B), meteorological forcing with atmospheric parameters varying in time (approach C).

Not considering the exchange process at the free surface can produce a good performance in the first layers of the reservoir (approach A). However, it may lead to overestimating the measured temperature in the deepest layers due to non-observed mixing processes. Moreover, it does not allow to simulate daily temperature oscillations.

Approach B with constant values of air temperature and wind, is the one that best represents the thermal stratification in terms of thermal distribution in the reservoir; however, there are excessive heat losses, which do not occur in reality. To better control these losses, Mesquita et al. (2020) recommend explicitly simulating evaporation in the reservoir.

Although approach C is the only approximation that shows hourly fluctuations in the temperature of the first layer of the reservoir, these fluctuations are larger than the observed ones. Simulation results show oscillations between 2 and 4.5 °C, whereas observed variations range between 1 and 2 °C. This may be due to the fact that the reservoir is located in a relatively sheltered canyon, as mentioned by Plec et al. (2021).

The results show that the parameter that mainly influences the thermal exchange model, thus improving the numerical reproduction of the thermal stratification processes, is the exchange coefficient A. A represent processes such as radiation, air convection in contact with water, and latent heat. Its correct determination improves the simulation results. According to the results of this study, when taking into account the exchange process between the water and the atmosphere under different forcings of the meteorological variables, a range of variation of A∈ [100, 900] W/m2/°C is best when the air temperature remains constant over time and lower values of A∈ [30, 40] W/m2/°C are best when weather conditions vary over time.

Although some uncertainties are present in the current numerical study about the A exchange coefficient tested in the simulations, the results are in a satisfactory agreement to data, especially considering the complexity of the studied phenomenon, as well as the equally complex dynamics of the various processes' interactions involved in its numerical simulation. In particular, regarding the water temperature at different depths, the following convergence results are observed: approach A maximum deviation from data reaching 3 °C, approach B maximum deviation from data reaching 4 °C, and approach C maximum deviation from data reaching 8 °C.

In conclusion, it can be stated that if site conditions and air–water temperature differences are known, A coefficient can be used as an approximate method of estimating the rate of exchange of natural water surfaces in a reservoir on an hourly basis during every period in the year.

Overall, we showed the potential and limitations of present three-dimensional hydraulic models in predicting the vertical distribution of temperature in tropical reservoirs. Future works will require rethinking the exchange formulae for the tropical context including evaporation, as the existing ones are not able to reproduce the variations in the water column of the reservoirs.

The authors would like to thank the Institut de Mécanique des Fluides de Toulouse, especially Maxime Pigou from the code and numerical simulations service and University of Medellin for providing us with the resources and support we needed to complete this project. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 22AP22035.

J.A. was involved in conceptualization, methodology, software, investigation, writing – original draft. H.R. was involved in validation, resources, writing – review & editing, supervision. L.C. was involved in validation, resources, writing – review & editing, supervision. T.B. was involved in validation, resources, writing – review & editing, supervision. J.E. was involved in conceptualization and supervision. L.J.M. was involved in data curation and supervision.

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

The authors declare there is no conflict.