Vertical Boundary Mixing Events during Stratification Govern Heat and Nutrient Dynamics in a Windy Tropical Reservoir Lake with Important Water-Level Fluctuations: A Long-Term (2001–2021) Study

Abstract

Physical processes play important roles in controlling eutrophication and oligotrophication. In stratified lakes, internal waves can cause vertical transport of heat and nutrients without breaking the stratification, through boundary mixing events. Such is the case in tropical Valle de Bravo (VB) reservoir lake, where strong diurnal winds drive internal waves, boundary mixing, and hypolimnetic warming during stratification periods. We monitored VB during 21 years (2001–2021) when important water-level fluctuations occurred, affecting mixing and nutrient flux. Stability also varied as a function of water level. Hypolimnetic warming (0.009–0.028 °C day⁻¹) occurred in all the stratifications monitored. We analyzed temperature distributions and modeled the hypolimnion heat budget to assess vertical mixing between layers (0.639–3.515 × 10⁻⁶ m³ day⁻¹), vertical diffusivity coefficient K_Z (2.5 × 10⁻⁶–13.6 × 10⁻⁶ m² s⁻¹), and vertical nutrient transport to the epilimnion. Nutrient flux from the metalimnion to the epilimnion ranged 0.42–5.99 mg P m⁻²day⁻¹ for soluble reactive phosphorus (SRP) and 5.8–101.7 mg N m⁻²day⁻¹ for dissolved inorganic nitrogen (DIN). Vertical mixing and the associated nutrient fluxes increase evidently as the water level decreases 8 m below capacity, and they can increase up to fivefold if the water level drops over 12 m. The observed changes related to water level affect nutrient recycling, ecosystemic metabolic balance, and planktonic composition of VB.

1. Introduction

Understanding mechanisms that control eutrophication and oligotrophication is an issue that involves fascinating theoretical and applied possibilities. The vertical exchange between the warm surface layer (epilimnion) and the colder hypolimnion below the thermocline boundary in stratified systems is an understudied critical process that may affect the nutrient dynamics and, consequently, their trophic condition [1]. In lakes characterized by stable stratification, vertical transport of dissolved materials is generally restricted by the metalimnetic barrier. Therefore, in stratified eutrophic lakes, primary producers inhabiting the epilimnion are often limited by the low availability of dissolved nutrients, while in the hypolimnion, concentrations of these elements can be much higher due to the microbial transformation of settling organic material and to their release from the bottom sediments, e.g., [2]. Thus, any physical mechanism responsible for significant cross-metalimnetic upward mass transport will be important in controlling the system’s productivity and its trophic condition.

External loading and deep-water entrainment are usually the two dominant sources of nutrients to lakes at an annual time scale, yet their relative importance during the period of stratification and low nutrient availability is under discussion [3]. As crossing the metalimnetic density barrier requires energy, entrainment—the transport of nutrients from the deep water—has been considered low or negligible during strong stratification [4]. Recent works, however, suggest that the vertical flux of nutrients across the thermocline may be an important driver of epilimnetic metabolism in some stratified lakes [5]. In this case, the energy required has often been shown to come from processes such as wind mixing and shear production within boundary layers of the lake [6]. Intense periodical winds can drive internal waves, which can cause increased dissipation of turbulent kinetic energy, both at the periphery and the center of the lake [7]. A recent study [8] reports that the synchronization of internal waves can even provoke enhanced shear in the hypolimnion and vertical diffusivity (Kz), reaching up to 10⁻⁴ m² s⁻¹. Nevertheless, processes driving the enhancement of boundary mixing events and hypolimnetic entrainment are still understudied, and even more seldom approached in tropical lakes [9,10].

This kind of process likely occurs in Valle de Bravo (VB) [11], a deep monomictic tropical reservoir lake in Central Mexico, which is daily swept by strong diurnal (12:00 a.m.–19:00 p.m.) winds (7.4–16.5 m s⁻¹). Important vertical diurnal displacements of the thermocline at the extremes of the reservoir lake reveal the internal waves driven by diel winds in VB [12]. Through the interaction of the thermocline with the bottom and margins of the lake, these internal waves drive daily boundary mixing events in VB during the stratification period, when hypolimnetic warming occurs [11]. Regular, long-term monitoring of VB has shown that this hypolimnetic warming occurs during every seasonal stratification period. While boundary-mixing events are not strong enough to break the stratification in VB, their cumulative effect drives this observed increase in the hypolimnetic temperature, which varies depending on the intensity and frequency of the mixing events at VB [11,13,14].

Additionally, because VB is a reservoir used for water supply, it exhibits important water level fluctuations (WLFs), as a result of the variations on water demand and on rainfall in its basin [15]. This is a relevant characteristic of VB because it is well known that important WLFs can affect multiple processes in water bodies [9,13,16]. In particular, long-term studies on the metabolic balance of VB have shown that it varies as a function of water level and suggest that low water levels probably intensify the boundary mixing events that occur in VB [14].

In fact, it is well known that mixing processes and their variability can significantly affect multiple aspects, including nutrient availability and budgets [5,7,9,11,15], ecosystem metabolism [5,14], as well as plankton composition and dominance [13,17] in all kinds of water bodies. In VB, there are, additionally, other important management issues—such as toxic cyanobacterial blooms—that challenge its value and use as a touristic hot spot [11], which are also affected by these mixing processes and their variability associated with WLFs. This has been outlined as a broadly found issue by [9], who stressed that the response of aquatic ecosystems, particularly deep lakes, to WLFs is an understudied field of crucial importance to the management of water resources.

As important WLFs have been occurring in VB in the past two decades, their long-term study can contribute to the understanding of their effects on mixing and biogeochemical processes. Such an understanding can be useful to cope with the eutrophication of deep monomictic lakes and reservoirs, which is an urgent management issue [18]. We monitored this reservoir lake for over 21 years, and in this study, we (1) report on the variability of temperature distributions in the water column of the windy reservoir lake of VB and (2) used the hypolimnetic temperature increase rates observed during each of the stratification periods to calculate the vertical mixing, the coefficient of vertical eddy diffusivity (Kz), and the vertical nutrient entrainment into the epilimnetic zone for each stratification of 2001–2021. As nutrients are non-conservative and rapidly affected by multiple processes (e.g., [19,20,21,22]), we did not compare these fluxes with the epilimnetic nutrient concentrations but rather with the changes in its biogeochemical processes and communities. Such a long-term series reveals the variation in the intensity of water mixing and nutrient entrainment as a function of water level fluctuations (WLFs). Comprehension of reservoir hydrodynamics is essential for maintaining good water quality and drinking water safety, as well as for predicting future changes [23]. Furthermore, the important ecological effects and implications of vertical mixing variations are discussed so that they can be used, along with other relevant considerations, in the context of the management of tropical lakes and reservoirs [18].

2. Materials and Methods

2.1. Study Area

Valle de Bravo (VB) has a surface of 18.6 km², a mean depth of 21.1 m, and a maximum depth of 38.6 m. It is the largest of seven reservoirs in central Mexico (19°21′30″ N; 100°11′00″ W) that form the Cutzamala System, which provides 30% of the drinking water supply to the Mexico City Metropolitan Area and other cities. The maximum capacity of this reservoir is 391 × 10⁶ m³. There are no spill outflows, but water is regularly withdrawn from VB by the National Water Agency (CONAGUA); on average, 221 × 10⁶ m³ have withdrawn annually from this reservoir. Water is replenished to this high-altitude (1830 m a.s.l. at its maximum water level) tropical reservoir through its 547 km² watershed by four main rivers (Amanalco, Molino, González y Carrizal) and at least three sewage outlets, including Tizates, a former river [15].

Through these rivers and sewages, VB receives very high external loadings of P and N (120.8 × 10³ kg P y⁻¹ and 591.8 × 10³ kg N y⁻¹ during 2002–2005, as already reported [15]). Due to these high external loads, this reservoir lake has become a hypereutrophic system, in which there are very high primary production rates (3.60 g C m⁻² day⁻¹ measured [14] during 2006–2015) in the photic layer, at the same time when the reservoir lake as a whole is heterotrophic. Due to its depth, VB behaves as a monomictic lake [11], and when it stratifies, it exhibits drastic biochemical differences between the epilimnion and the hypolimnion. During stratification, the hypolimnion becomes anoxic, and nutrient concentrations in the hypolimnion are high (SRP = 1.25 $\mathsf{\mu }\mathrm{M}$, DIN = 36.3 $\mathsf{\mu }\mathrm{M}$, see [11]), while in the epilimnion, they decrease to an order of magnitude lower levels in the epilimnion (SRP = 0.28 $\mathsf{\mu }\mathrm{M}$, DIN = 3.1 $\mathsf{\mu }\mathrm{M}$, see [11]).

VB is located in the tropical latitudinal zone with trade winds, rainy around summer, and mostly dry the rest of the year. Climate plays an important role for VB, mainly due to rain and wind patterns. It is sub-humid, warm to temperate with very pronounced dry and rainy seasons. Mean annual precipitation is 947 mm, most of which (809 mm, 86%) fall in the rainy months (June–October), while barely 14% (138 mm) of the rain falls distributed throughout the six months (November–May) of the dry season. The mean temperature in VB is 18.7 °C, the minimum annual temperatures (mean 9.0 °C) occur during January, and the maximum (mean 28.0 °C) during May, before the onset of the rainy season. The diurnal temperature variation is important (mean 11.5 °C) and can reach up to 14 °C in January. The climate of this deep reservoir lake is also characterized, as previously mentioned, by strong diurnal winds that blow along its two arms, of which the longer has a 6.9 km length (Figure 1). The wind at VB exhibits a regular diurnal pattern (12:00 a.m.–19:00 p.m.) blowing from the reservoir’s dam (WNW) into the valley arms (SSE). Diurnal wind speed averages 5.5 m s⁻¹ and frequently reaches over 10 m s⁻¹ [11].

Limnologically, VB exhibits a warm monomictic thermal regime, similar to other water bodies of the Mexican tropical highlands, except for the occurrence of a continuous increase in its hypolimnetic temperature during the stratification period [11]. In contrast, no hypolimnetic temperature increase has been observed or reported in any other water bodies of the Mexican tropical highlands, except for an almost negligible increase in [24].

2.2. Long-Term Monitoring of VB

VB monitoring was performed monthly at up to 17 sampling stations distributed throughout the water body and along the wind axis (Figure 1). During the initial years of the monitoring (2001–2003), all the stations shown in Figure 1 were sampled to identify and include any spatial variations that could occur. However, as the existence of horizontal homogeneity in VB was confirmed [11,25], the number of stations was reduced, initially to five stations, and from September 2007 on to a single sampling station at the center of the lake, to extend the viability of the long-term monitoring that has proceeded since then. The reservoir’s water level during each sampling date was obtained from the government official reports, and the area of the reservoir lake at different depths was obtained from the bathymetry reported by [11].

Temperature and dissolved oxygen (DO) were measured at 1 m vertical intervals using an Hydrolab DS4/SVR4 (Austin, TX, USA) (from February 2001 to September 2004) or a Yellow Springs Instruments (YSI, Yellow Springs, OH, USA) field probe (from October 2004 to January 2019) and an in situ Aqua Troll 500 (In-Situ, Fort Collins, CO, USA) (from February 2019 to September 2021), which were previously intercalibrated. Water samples were collected with a 1.5 L Uwitec (Mondsee, Austria) sampler at regular vertical intervals (0, 1, 2, 4, 8, 12, 17, 20, 24, 28, and 32 m, when possible) all the way down to the bottom of the lake. Subsamples for nutrients (dissolved inorganic nitrogen (DIN) = [NH₄⁺ + NO₂⁻ + NO₃⁻], and soluble reactive phosphorus (SRP) analysis were held in polypropylene containers after filtration with 0.45 and 0.22 µm Millipore™ type HA (Bedford, MA, USA) nitrocellulose membranes and were kept frozen until their analysis (within 24 to 48 h) in a San plus Skalar (Breda, The Netherlands) segmented-flow continuous auto-analyzer, using the methods and circuits suggested by [26].

2.3. Vertical Mixing Model and Calculations

Our approach is based on the increasing evidence [6,27,28,29] that turbulence generated by internal waves breaking in boundary mixing events can transport nutrients and gases vertically through the metalimnion in systems where a source of external energy, such as wind, drives these processes, as summarized in Figure 2. In VB, these mixing events drive continuous warming of the hypolimnion that continues along all the stratification periods.

To assess if hypolimnetic warming reduced the stability of the stratification (S), stability was calculated for each sampling date following [30], integrating from the surface to the maximum depth (Equations (1) and (2)). Stability is expressed in work units per surface lake area because it represents the amount of work that would be required to mix the reservoir lake into an isothermal state. It describes how much the center of mass of the lake has been lowered by the stratification process.

S = g/A_o ∑ (ρ_z − ρ_m) (z − z_g) A_zΔz

(1)

where

ρ_m = 1/V ∫A_zρ_zdz

(2)

where S is the stability of stratification (J m⁻²); g is the acceleration of gravity (m s⁻²); A_o is the surface area of the reservoir lake (m²); ρ_z is the water density at depth z (g m⁻³); ρ_m is the mean water density of the lake (g m⁻³); z_g is the depth at which the mean density is found (m); A_z is the area at any depth z; V is the volume of the reservoir lake (m³).

2.4. Hypolimnetic Heat Balance and Calculation of Vertical Water Mixing

There are not any geothermal or groundwater sources of heat in VB that could be driving the warming of the hypolimnion [11]. The high turbidity of the reservoir water—mean Secchi depth is of only 1.4 m (stratification averages range: 1.15–1.65) during the stratification periods [14]—constrains the penetration of solar radiation heating to the first few meters of the epilimnion. Therefore, the hypolimnetic thermal increase we found was entirely attributed to vertical mixing events. We used the simplest model possible, considering the three distinct water layers in a monomictic lake: the warm epilimnion, metalimnion or thermocline layer, and the colder hypolimnion (Figure 2). As there are no other heat inputs to the hypolimnion of VB, it is reasonable to trace back the vertical mixing from the metalimnion to the hypolimnion, Mix_meta-hypo, from the observed hypolimnetic temperature increases. This was carried out through a heat balance of the hypolimnion (Equation (3)), considering both the heat exchange due to a) vertical mixing of water (ΔQhypo_mix, cal d⁻¹) and b) thermal conductivity flux (ΔQhypo_cond, cal day⁻¹), as

ΔQhypo_tot = ΔQhypo_mix + ΔQhypo_cond

(3)

Total heating of the hypolimnion, ΔQhypo_tot, was calculated (Equation (4)) for each period between samplings (Δt, days), using the observed temperature increase (ΔT_hypo, °C), the calorific capacity of water, (c, cal g⁻¹ °C⁻¹), the volume of the hypolimnion (V_hypo, m³), and the mean hypolimnetic density (ρ_hypo, g m⁻³), for a period Δt, as

ΔQhypo_tot = (ΔT_hypo c/Δt) × V_hypo × ρ_hypo

(4)

The conductivity heat flux, ΔQhypo_cond, was calculated from the vertical temperature gradient (ΔT_epi-hypo) at the metalimnion–hypolimnion frontier and the water thermal conductivity. Then, the mixing heat exchange (ΔQhypo_mix) was calculated by subtracting the thermal conductivity flux from the total heat exchange (resolving from Equation (3)).

In turn, the volume of water exchanged (Mix_meta-hypo, m³ d⁻¹) was calculated from the obtained heat exchange due to mixing (ΔQhypo_mix, cal day⁻¹), as

Mix_meta-hypo = ΔQhypo_mix/ρ_hypo c (T_epi − T_hypo)

(5)

where: c is the calorific capacity of water (cal g⁻¹ °C⁻¹), ρ_hypo is the mean hypolimnetic density (g m⁻³), calculated as the mean of the ρ values measured in the hypolimnion, T_meta is the mean metalimnetic temperature (°C), and T_hypo is the mean hypolimnetic temperature (°C), all three calculated averaging the measurements from the two samplings used for each period calculation.

The vertical diffusivity coefficient, Kz (m² s⁻¹) sensu [6,8] was calculated from the vertical water flow (Mix_meta-hypo), multiplying it by the mean vertical distance between the two layers (Z_meta-hypo ~4 m) and dividing it by the exchange area (Area _meta-hypo) at the depth of the base of the metalimnion and the top of the hypolimnion (Equation (6)), which was estimated at 1.2 × 10⁷ m² using the bathymetry of the reservoir published in [11] and the corresponding conversion of units.

Kz = Mix_meta-hypo (Z_meta-hypo)/Area _meta-hypo

(6)

2.5. Calculation of Vertical Nutrient Fluxes

To assess the vertical flux of nutrients toward the epilimnion, estimates of both the metalimnetic concentration of nutrients and of the vertical rate of metalimnetic water mixed into the epilimnion (Mix_meta-epi) are needed.

Despite the hypolimnetic warming, because the thermal structure did not change notably in VB between samplings, and the reservoir lake remained stratified during the spring to autumn months, it is reasonable to assume a steady-state condition of vertical mixing, in which vertical water exchange between the epilimnion and the metalimnion is similar to vertical water exchange between this thermocline layer and the hypolimnion in both directions (Figure 2, Equation (7)).

Mix_epi-meta ≈ Mix_meta-hypo ≈ Mix_hypo-meta ≈ Mix_meta-epi

(7)

We, therefore, used the calculated Mix_meta-hypo to estimate Mix_meta-epi, the vertical water flux responsible for transporting nutrients to the epilimnion. Therefore, the nutrient entrainment, or flux to the epilimnion (F (nutrients)), for each period between samplings was calculated as the product of mean nutrient concentration in each layer multiplied by Mix_meta-hypo.

F (nutrients) = Mix_meta-hypo × (nutrients)

(8)

where F (nutrients) is the ascending nutrient flux (kg day⁻¹), Mix_meta-hypo is the exchange of water between layers (m³ day⁻¹), and [nutrients] is the mean metalimnetic concentration of SRP, ammonia, or DIN (kg m⁻³) (Equation (8)).

2.6. Data Statistical Analysis

Regression analysis and Pearson´s determination coefficients were used to investigate the relations among the variables studied and to obtain the mean hypolimnetic warming rate for each of the stratification periods. Microsoft Excel for Mac 2019 (King County, WA, USA) was used for these analyses.

3. Results

3.1. Water Level Fluctuations (WLFs) in VB during 2001–2021

A seasonal variation on the water level of VB was observed during the full study period (Figure 3a), which followed the seasonal rain variations for the region, where the rainy season is from June to October. These seasonal water level variations were generally of only 4–6 m, between the reservoir lake’s capacity (maximum level 1830 m a.s.l.) at the end of the rainy season and a minimum of 1824 to 1826 m a.s.l. by the end of the dry season (Figure 3a). However, the level of the reservoir lake varied more drastically during 2005–2013, decreasing over 12 m below its capacity and remaining well below its maximum level for several years, except for 2010 (Figure 3a). Then, VB returned to water levels close to its capacity for several years (2014–2020), until the level fell drastically again in 2021. Overall, four extraordinary minimum water levels were reached in the summers of 2006 (1819.7 m a.s.l.), 2009 (1817.9 m a.s.l.), 2013 (1818.3 m a.s.l.), and 2021 (1816.2 m a.s.l.), which are outlined in red color in Figure 3 and Figure 4 to facilitate their identification. Despite these high WLFs, stability maintained marked annual cycles that show that the stratification was not broken during any of the stratification periods of 2001–2021 (Figure 3b).

3.2. Thermal Regime and Stratification

VB behaved as a warm monomictic lake every year during 2001–2021, with a mean stratification period of 219.8 ± 30.6 days. A typical seasonal pattern occurred in the reservoir during the 21 years studied: the maximum surface temperature was reached in June, with slight variations (23.0–23.9 °C) among years (Figure 4a). Circulation began in November, and the water column remained homogeneous every winter, at least from December through January. The time variations of the vertical thermal distribution show that the epilimnion extended quite regularly among years to a depth of ~6 m, and a little deeper in some of the years, particularly in the two lowest level years (2009 and 2021). As expected, the stability of the water column also followed a seasonal pattern, reaching high values (>20.0 J m⁻²) during all the stratification periods and decreasing to zero every year during the mixing periods (Figure 3b). This seasonal pattern of stability confirms that stratification was never broken during any of the stratifications between 2001 and 2021. Nevertheless, the mean stability observed during each of the stratifications diminished significantly as the water level decreased (Figure 3c). In fact, as shown in Figure 3c, stability exhibited a direct and significant polynomial relationship with the water level of the reservoir lake.

3.3. Hypolimnetic Warming

The mean temperature of the hypolimnion of VB increased between consecutive samplings during each and every stratification period of the 21 years here reported (Figure 4). The rates of temperature increase between consecutive samplings ranged overall 0.0003–0.04 °C day⁻¹, while the mean rates for each entire stratification period, as estimated by linear regressions in Figure 4b ranged from 0.009 to 0.028 °C day⁻¹, and the maximum mean increases occurred during the stratifications of 2021 (0.028 °C day⁻¹), 2009 (0.022 °C day⁻¹), 2013 (0.020 °C day⁻¹), and 2006 (0.019 °C day⁻¹). Although linear regression yielded high coefficients of determination for all the stratification periods (r² ranging from 0.936 to 0.999), some of the stratifications showed a slight sigmoidal trend (e.g., 2016, 2012, 2018, 2005), which could be related to the transitions between stratification and circulation periods. The mean thermal increase during the stratification showed a strong (r² = 0.897) second-order inverse relation with the mean water level for that year: the lower the level of the lake, the higher the rate of hypolimnetic warming (Figure 4c).

3.4. Mix_meta-hypo, Kz, and Nutrient Fluxes

The rate of vertical water mixing (Mix_meta-hypo) ranged from 0.639 × 10⁶ m³ day⁻¹ to 3.515 10⁶ m³ day⁻¹ among the stratifications, with lowest in 2015 and highest in 2021 (

Table 1. Mean vertical water mixing rates (10⁶ m³ day⁻¹), mean vertical diffusivity coefficients, Kz (m² s⁻¹), and mean rates of nutrient vertical flux, F, of SRP (mg P m⁻² day⁻¹) and DIN (mg N m⁻² day⁻¹) calculated for the stratification periods of 2001–2021 in VB.

Stratification Period	Mixmeta-hypo (106 m3 day−1)	Kz (m2 s−1)	F (SRP) (mg P m−2 day−1)	F (DIN) (mg N m−2 day−1)
2001	0.861	3.3 × 10−6	0.79	15.9
2002	0.924	3.6 × 10−6	2.03	28.0
2003	1.057	4.1 × 10−6	0.57	20.3
2004	0.880	3.4 × 10−6	1.03	21.7
2005	1.046	4.0 × 10−6	0.68	21.1
2006	1.691	6.5 × 10−6	1.12	71.5
2007	1.212	4.7 × 10−6	0.42	13.0
2008	1.963	7.6 × 10−6	4.75	31.6
2009	3.373	13.0 × 10−6	5.99	97.1
2010	1.439	5.6 × 10−6	2.64	17.0
2011	1.271	4.9 × 10−6	1.54	18.6
2012	1.446	6.6 × 10−6	1.95	30.4
2013	2.349	9.1 × 10−6	5.11	73.8
2014	0.832	3.2 × 10−6	2.24	7.6
2015	0.639	2.5 × 10−6	1.30	5.8
2016	0.712	2.7 × 10−6	1.52	12.5
2017	1.296	5.0 × 10−6	1.17	32.0
2018	1.180	4.6 × 10−6	0.56	30.4
2019	1.059	4.1 × 10−6	2.65	36.3
2020	0.691	2.6 × 10−6	1.20	41.0
2021	3.515	13.6 × 10−6	5.27	101.7
2001–2021 Mean	1.289	3.1 × 10−6	2.12	34.6

). These mixing rates correspond to a daily volume of 0.16% to 0.90% of the 391 × 10⁶ m³ day⁻¹ total capacity of the reservoir lake (see Section 2.1. Study Area). Furthermore, these rates showed an inverse relationship with the mean water level for each year, with a good fit (r² = 0.893) through a polynomial second order regression (Figure 5a).

The mean calculated vertical diffusivity coefficients (Kz) followed a similar pattern among years, ranging from a 2.5 × 10⁻⁶ m² s⁻¹ average during the stratification of 2015 to a mean of 13.6 × 10⁻⁶ m² s⁻¹ during the stratification of 2021 (Table 1, Figure 5b) when the water level was lowest during the full period. On the monthly scale of our calculations, the range of the Kz values calculated was wider, from as low as 4.0 × 10⁻⁸ m² s⁻¹ to as high as 5.0 × 10⁻⁵ m² s⁻¹ (

Table 2. Mean and range of the vertical diffusion coefficients Kz (m² s⁻¹) obtained for VB during the stratification periods of 2001–2021 and a comparison with Kz reported for other water bodies. The mean and range (in parenthesis) of Kz values are included. Data are ordered from higher to lower Kz mean values, and the position of VB among other lakes and reservoirs is outlined in bold characters and green shading. When not reported, mean Kz were calculated from the values available. The value for molecular conductivity is included to show the lower limit of Kz when mixing is null, and only molecular movements occur.

Kz		Location	Reference
Mean	(Range)
9.9 × 10−6	(1.2 × 10−7–5.0 × 10−5)	VB (monthly)	This study
3.1 × 10−6	(2.5 × 10−6–13.6 × 10−6)	VB (annual mean)	This study
5.0 × 10−4	(1.0 × 10−5–1.0 × 10−3)	North Lake, Australia	[31] MacIntyre 1993
1.0 × 10−4	(8.0 × 10−5–1.3 × 10−4)	Lake Rotowhero, New Zealand	[32] Brookes et al. 2013
6.5 × 10−5	(2.2 × 10−5–1.1 × 10−4)	Lake Lanao, Philippines	[33] Lewis 1982
5.1 × 10−5	(1.0 × 10−6–1.0 × 10−4)	Lake Biwa, Japan	[8] Auger et al. 2014
5.1 × 10−5	(1.0 × 10−6–1.0 × 10−4)	Kranji Reservoir, Singapore	[34] Yang et al. 2019
4.5 × 10−5	(1.0 × 10−6–9.0 × 10−5)	Green Lake,	[35] Stauffer 1987
4.5 × 10−5	(8.0 × 10−5–1.0 × 10−5)	Lake Kinneret, Israel	[36] Saggio and Imberger 2001
1.8 × 10−5	(1.6 × 10−5–2.1 × 10−5)	Fitzroy River Barrage, Australia	[37] Bormans et al. 2004
1.2 × 10−5		Lake Constance, Germany-Switzerland	[38] Preusse et al. 2010
9.9 × 10−6	(4.0 × 10−8–5.0 × 10−5)	Valle de Bravo (monthly)	This study
6.0 × 10−6	(9.0 × 10−6–1.1 × 10−5)	Lake Valencia, Venezuela	[39] Lewis 1983
6.0 × 10−6	(6.0 × 10−6–5.0 × 10−5)	Lake McIlwaine, Rhodesia	[40] Robarts and Ward 1978
5.0 × 10−6		Lake Matano, Indonesia	[41] Katsev et al. 2010
3.4 × 10−6	(1.0 × 10−7–2.6 × 10−5)	Mono Lake, California	[42] Jellison and Melack 1993
3.1 × 10−6	(2.5 × 10−6–13.6 × 10−6)	Valle de Bravo (annual mean)	This study
2.2 × 10−6	(3.0 × 10−6–8.0 × 10−5)	Lake Alpnach (central Switzerland)	[43] Goudsmith et al. 1997
1.3 × 10−6	(2.0 × 10−7–2.5 × 10−6)	Castle Lake, California	[44] Jassby and Powell 1975
1.0 × 10−6	(1.0 × 10−8–1.0 × 10−4)	Mono Lake, California	[6] MacIntyre et al. 1999
1.0 × 10−8		Molecular conductivity	[6] MacIntyre et al. 1999

).

Vertical nutrient fluxes toward the epilimnion (Table 1) also varied inversely with the lake´s water level, following a pattern similar to that of Mix_meta-hypo and Kz (Figure 6). SRP vertical flux to the epilimnion ranged from 0.42 to 5.99 mg P m⁻² day⁻¹ among the different stratification periods and averaged 2.12 mg P m⁻² day⁻¹ for the full 21 years studied. DIN vertical flux to the epilimnion ranged from 5.79 to 101.7 mg N m⁻² day⁻¹ among the different stratifications and averaged 34.63 mg N m⁻² day⁻¹ for the 21 years studied. Significant inverse correlations (r² = 0.648 for SRP and r² = 0.845 for DIN) were obtained with the water level, although—as expected—nutrient fluxes exhibited a higher variability because of the additional variability of nutrient concentration, as can be seen in Figure 6.

4. Discussion

The access of nutrients into the photic zone has major implications for the primary production and trophic webs of water bodies. This is largely controlled by physical factors that include the wind regime and water level fluctuations. The approach used here was derived from long-term observation throughout two decades of monthly monitoring at VB, where strong winds, eutrophication, and water extraction for Mexico City interact in a system where explosive phytoplankton blooms are frequent [11,17]. In such a complex system, the assessment of the vertical exchange of water and nutrients and its variations can help to reveal the mixing and nutrient environment to which plankton is exposed, and nutrient budgets can be refined. From this understanding, insight and prediction of the processes involved could be derived, which, eventually, become useful tools for the management of similar lakes and reservoirs in the global change scenarios to come.

Even though major WLFs have been recorded in VB during the past 21 years, our observations demonstrated that this tropical reservoir lake behaves regularly as a monomictic tropical lake, exhibiting a long and stable stratification from spring to autumn and a circulation period during the winter [10,11,14,15]. Moderate seasonal changes in the water level of VB occur every year as a function of seasonal rainfall variations in its basin, but higher WLFs are results of longer-term climatic variability and management decisions that consider requirements of water for human consumption and water availability within the Cutzamala water supply system as a whole; thus, they occur on multi-annual time scales. Therefore, a long-term series, as the one here reported, is required to include enough water level variability to allow the identification of patterns and the effects of these variations in tropical deep lakes such as VB. The magnitude of WLFs in VB during this period was so high that a level decrease of up to almost 14 m was registered, which is more than half of the mean depth of the reservoir lake when it is at its maximum level. It is, however, remarkable that VB remained stratified during the warm months (March to November) of each and all of the years (2001–2021) included in our sampling period.

4.1. Hypolimnetic Warming Allows the Assessment of Vertical Mixing Rates during Stratification

Our temperature measurements also showed that the hypolimnion warmed steadily during every stratification in VB throughout the entire observation period, although the warming rate varied from year to year. We attribute this warming to vertical mixing events that are likely driven by internal waves caused by the wind pattern in VB, which has a very regular diurnal oscillation [11]. This is in agreement with the findings of [6,27,28,29] in similar systems, where internal waves have been proposed to drive boundary mixing and generate vertical nutrient transport. In the absence of other heat sources than vertical mixing in VB, we used the hypolimnetic warming rates to assess vertical mixing in VB during the stratification periods. Table 2 shows that our calculated coefficients (mean 9.9 × 10⁻⁶ m² s⁻¹, range 1.2 × 10⁻⁷ m² s⁻¹ to 5.0 × 10⁻⁵ m² s⁻¹) are amidst the range of values reported in specific studies of vertical diffusivity, even on the lower part of the range but always above molecular conductivity (10⁻⁸ m² s⁻¹), supporting the fact that we were investigating mixing processes. Higher average Kz than those we calculated have been reported, both in tropical (e.g., Lake Lanao [33] in the Philippines and Kranji Reservoir [34] in Singapore) and in temperate (e.g., Lake Constance [38] in the border between Germany and Switzerland) systems. In particular, at the intensely studied Lake Kinneret, where internal waves and boundary mixing are known to drive vertical mixing [36], measured vertical diffusivities are on average about an order of magnitude higher than our mean of monthly Kz estimates.

Nevertheless, our monthly vertical diffusivity values are, on average, comparable and slightly higher than those reported for other tropical lakes, including Lake Valencia [39] in Venezuela, Lake McIlvaine [40] in Rhodesia, and Lake Matano [41] in Indonesia, and are also higher than those found in some temperate lakes, such as Mono Lake [6,42], Lake Alpnach [43] in Switzerland and Castle Lake [44] in California. Similar to our findings, Kz values have been found to vary around three orders of magnitude within the metalimnetic layer [45].

When comparing with other vertical diffusivity reports, it should be noted that our estimates assessed the combined effect of vertical mixing events that occurred during a full month, i.e., between two consecutive sampling dates, while direct diffusivity measurements (such as most of those pooled together in Table 2) are performed on much smaller time scales and likely register the peak values that occur during mixing events. In contrast, our estimates are smoothed values over the full period between sampling and are likely more representative of the long-term mixing rates. For instance, at Lake Kinneret, where mixing processes are amongst the most well known, a study [36] outlines the variability that might affect the representativity of short-term measurements: “most of the time the vertical flux through the metalimnion was negligible, but, at times, the eddy diffusivity did reach values as high as 10⁻² m² s⁻¹”. We believe that, by integrating the effect of these relevant vertical mixing pulses, our results likely better represent actual mixing rates at the near-monthly scale. Nevertheless, to fully verify this in future research, we will explore the mixing intensity at short-time scales.

4.2. Vertical Mixing Increases Up to Fivefold as the Water Level Decreases

In terms of their importance for VB, the magnitude of the vertical water exchange rates here calculated (Table 1) imply that 5–11% of the total volume of the reservoir lake is being exchanged monthly between layers during the stratification period of the years when the water level is relatively high. However, as Figure 5 shows, vertical mixing varied significantly among the stratification periods as a function of the lowest water level for each year. During the lowest water level years, the vertical exchange can reach up to 26% (e.g., 2009 and 2021) of the lake’s volume in a single month. When integrated for the full stratification periods, these volumes of water represent 40–83% of the total lake´s volume being exchanged between layers, and this exchange can reach up to 190% of the lake´s volume in the lower level and stronger mixing years, such as 2009 and 2021.

These magnitudes suggest that vertical transport can be much higher than what would be expected from the traditional paradigm of stratification and metalimnetic barrier [2] and could also provide a physical mechanism to explain the sustained heterotrophic behavior (e.g., sustained high respiration rates) found for VB even under stratification conditions [14], as discussed further ahead.

The calculated vertical nutrient fluxes also increased substantially during the low water level stratifications in VB (Figure 6), reaching the highest values among the relatively few systems for which these fluxes were assessed (

Table 3. Mean Nutrient fluxes obtained for VB during the stratification periods of 2001–2021 and during the stratifications with the lowest water level. Fluxes reported from other lakes are also included for comparison. Other lakes’ data are ordered by increasing SRP flux, and the position of VB fluxes among these values is outlined in bold characters and green shading.

SRP (mg P m−2 day−1)	Ammonia (mg N m−2 day−1)	DIN (mg N m−2 day−1)	Lake	Reference
0.10			Lake Fish, Wisconsin, USA	[35] Stauffer 1987
1.60	6.2		Lake McIlwaine, Rhodesia	[40] Robarts and Ward 1978
2.12	13.9	34.6	Valle de Bravo (overall mean)	This study
2.71			Lake Kivu, Africa	[46] Pasche et al. 2009
2.59			Lake Mendota, Wisconsin, USA	[47] Kamarainen et al. 2009
~3.00	~20.0		Lake Malawi/Nyasa, Africa	[48] Hamblin et al. 2003
3.20			Lakes of Taylor Valley, Antarctica	[51] Edwards and Priscu 1995
5.50			Lake Green, Wisconsin, USA	[35] Stauffer 1987
5.99	66.0	101.7	Valle de Bravo (lowest water level)	This study
15.20			Lake Delavan, Wisconsin, USA	[35] Stauffer 1987
	39.2		Mono Lake, California, USA	[6] MacIntyre et al. 1999

). When the average for all the years in the 2001–2021 period are considered, including both lower and higher-level years, our estimates of vertical nutrient fluxes of SRP and ammonia in VB are similar to those that have been calculated for other lakes (Table 3), such as Kivu [46], Mendota [47], McIlwaine [40] and Malawi [48]. However, the vertical nutrient fluxes for the stratification of the years with the lowest water level (2009 and 2021) would be higher in VB than all those in Table 3, with the only exception of Lake Delavan in Wisconsin, where an SRP flux as high as 15.2 of mg P m⁻² d⁻¹ was estimated by [35]. This can be due to, in part, the relatively low proportion of our estimated SRP flux with respect to the DIN flux (overall mean 36.2), which can be due to P retention in the sediments] [49,50], as previously suggested by [15]. Nevertheless, when low-water years are compared with the rest, our results clearly show that the vertical nutrient flux from the hypolimnion increases in VB by an order of magnitude when the water level drops sharply (≈12 m), a factor that has not been considered in water reservoir management until now.

4.3. The Contribution of Vertical Nutrient Fluxes to the Surface Layer Budget Can More than Duplicate as the Water Level Decreases

It is also interesting to compare the vertical transport here calculated (due to internal loading) with the high N and P external loadings VB receives [15], mainly from urban wastewaters and agricultural land uses in the VB basin [52]. The nutrient vertical entrainment to the euphotic zone here reported represents, on average, about 12% of the total phosphorus annually entering the reservoir through its main tributaries (16.9 ± 8.7 mg P m⁻² day⁻¹ [15]). Nevertheless, this proportion clearly increases when the level falls below 8 m (Figure 6), and it can more than duplicate, reaching up to >35% if the level decreases to 12 m or more.

Thus, the water level can alter significantly the amount of nutrients available in the euphotic zone of the reservoir lake during stratification, and therefore, its primary production and functional eutrophication level, as identified by [14]. This is consistent with the particularity of VB being more productive during stratification than during the circulation period [13,14,17], during which nutrients would not be segregated to the hypolimnion and circulate freely. In fact, VB reaches its highest phytoplanktonic diatoms biomass peaks during the stratification periods, contrasting with warm monomictic lakes and reservoirs in either temperate zones or high altitudes at tropical latitudes which typically reach their highest phytoplanktonic diatoms biomasses during circulation [17].

The intensification of vertical mixing and nutrient transport during low water stratifications could favor more intensive recycling of nutrients in VB when the water level is low, while N and P fluxes to the sediments through sedimentation would likely be higher when the water level is high, and vertical transport decreases. Furthermore, the changes in metabolism found in VB [13], associated with the variation in vertical mixing, have been found to alter both the metabolic balance and the role of the reservoir lake as a source or sink of carbon to the atmosphere [14], and all this would depend largely on the variation in water level.

4.4. The Assessment of Vertical Nutrient Fluxes Also Helps Understanding Some Ecological Variations

A cascade of ecosystemic processes can be affected by boundary mixing and vertical nutrient transport, as has been addressed broadly in the literature [53,54,55,56,57,58]. At the same time, the temperature and nutrient signals of these processes are rapidly blurred in the surface layer, both by surface heat transport and by the multiple processes that make nutrients non-conservative, including competitive assimilation and transformation processes (e.g., [19,20,21,22]). Furthermore, N and P nutrients reaching the epilimnion are readily assimilated, keeping its surface concentration low, for example, in a typical stratification in VB [11]. However, there are also multiple indirect effects of vertical mixing and nutrient supply that have been observed in VB. Previous reports of oxygen dynamics in VB and plankton community composition and abundance [17,25] have been correlated with WLFs and boundary mixing, and comprise indirect evidence consistent with the nutrient transport results presented here. In more detail, the observation of ephemeral patches of oxygen below the photic zone (following the proposal of [59]), in the typically anoxic summer–autumn hypolimnion of VB [13], also supports the occurrence of frequent (≈daily) events of mixing during stratification. Moreover, low water periods, when our results indicate higher vertical mixing rates, were also correlated with higher respiration rates in VB in a 10 year-time series [14], pushing the metabolic balance to ecosystemic heterotrophy. This is likely because enhanced mixing makes possible the oxidation of reduced species accumulated in the hypolimnion, including ammonia [11] and organic matter [15].

Boundary mixing events in the summer stratified period have been also correlated with important impacts in the plankton communities, for example, with a diversification and biomass increase in the zooplankton community [13]. In particular, cladocerans and copepods have been observed to increase their abundance during low water level periods [13]. This is remarkable because rotifers are normally the most prevalent zooplankton group in VB [60]. Additionally, the enhancement of phytoplankton biomass has also been reported during a low water level stratification compared with the same period in a year of high water level conditions [17]. This is consistent with—and would be expected to occur as a result of—the intensification of nutrient inputs to the epilimnion through boundary mixing here predicted. Furthermore, important changes in phytoplankton groups were also observed in the same study, in which low water levels were characterized by increasing biomass of planktonic diatoms and a significant decrease in harmful cyanobacteria (filamentous Nostocales and coccoids, e.g., Microcystis). The effects of these changes in phytoplankton [17,58] and zooplankton [60] in the rest of the ecosystem are still to be studied. In particular, changes in zooplankton composition could favor fisheries availability [13]. Finally, there are also important effects of WLFs on littoral communities [9], which are still to be studied in VB.

Comprehension of lake and reservoir hydrodynamics is key for good water quality, for drinking water safety, as well as for predicting future changes [23]. Understanding the relationships between temperature and vertical transport is key to modeling emergent ecosystem functions in a warmer world [18]. In VB, vertical mixing and nutrient fluxes increases are evident when the level decreases below 8 m and can increase several-fold when the level falls down below 12 m. Altogether, the multiple consequences of WLFs in windy tropical lakes such as VB can alter both the ecology of the reservoir lake and its suitability for sustained water exploitation by the main users. The effects we identified here should be accounted for by water managers [18,23], who should outweigh them against the particular needs of water distribution and other reservoir lake uses when making decisions that may affect the water level of reservoirs and lakes.

5. Conclusions

The monomictic reservoir lake of VB exhibited hypolimnetic warming during all the stratification periods of 2001–2021. Vertical mixing calculated from this warming showed there is a significant vertical exchange in this reservoir lake during the stratification months, comparable to the vertical diffusivity found for other tropical and temperate water bodies where vertical mixing is important. Our results showed that vertical mixing during stratification increases in VB as the level of the reservoir lake decreases. We attribute this to the enhancement of boundary mixing events as internal waves intensify their interaction with the lake´s bottom and margins, as found for other lakes and reservoirs that have a similar wind forcing. Vertical diffusivity coefficients calculated for VB are within the range found in other water bodies, although they average the vertical mixing occurring over monthly-to-annual scales, while most direct measurements reported for other water bodies address peak mixing pulses.

The vertical nutrient transport in VB would be enough to overcome the stratification barrier and sustain the observed eutrophic conditions during stratification. This is consistent with the fact that VB is more productive during stratification than in the circulation period when nutrients are not segregated to the hypolimnion. Mixing and mean vertical nutrient entrainment to the epilimnion in VB is similar to in other studied systems but becomes much higher during years with low water levels. These increases become evident when the water level decreases more than 8 m below the reservoir capacity and can increase several-fold when the level drops down below 12 m. It is important that managers support long-term monitoring and also consider that water level also affects significantly the planktonic communities and the metabolic balance of the ecosystem, as well as the recycling of the external nutrient loading to the reservoir lake and the fraction accumulated at the bottom of VB.

Overall, our results provide insight into the physical controllers of the trophic status in windy tropical water bodies with high water-level fluctuations. It allows for the prediction of vertical mixing and nutrient entrainment as a function of the water level of lakes and reservoirs, providing an additional tool for water authorities and lake or reservoir users to improve their management strategies for VB and for similar water bodies.