Nitrogen and Phosphorous Retention in Tropical Eutrophic Reservoirs with Water Level Fluctuations: A Case Study Using Mass Balances on a Long-Term Series

Abstract

Nitrogen and phosphorous loading drives eutrophication of aquatic systems. Lakes and reservoirs are often effective N and P sinks, but the variability of their biogeochemical dynamics is still poorly documented, particularly in tropical systems. To contribute to the extending of information on tropical reservoirs and to increase the insight on the factors affecting N and P cycling in aquatic ecosystems, we here report on a long-term N and P mass balance (2003–2018) in Valle de Bravo, Mexico, which showed that this tropical eutrophic reservoir lake acts as a net sink of N (−41.7 g N m⁻² y⁻¹) and P (−2.7 g P m⁻² y⁻¹), mainly occurring through net sedimentation, equivalent to 181% and 68% of their respective loading (23.0 g N m⁻² y⁻¹ and 4.2 g P m⁻² y⁻¹). The N mass balance also showed that the Valle de Bravo reservoir has a high net N atmospheric influx (31.6 g N m⁻² y⁻¹), which was 1.3 times the external load and likely dominated by N₂ fixation. P flux was driven mainly by external load, while in the case of N, net fixation also contributed. During a period of high water level fluctuations, the net N atmospheric flux decreased by 50% compared to high level years. Our results outlining water regulation can be used as a useful management tool of water bodies, by decreasing anoxic conditions and net atmospheric fluxes, either through decreasing nitrogen fixation and/or promoting denitrification and other microbial processes that alleviate the N load. These findings also sustain the usefulness of long-term mass balances to assess biogeochemical dynamics and its variability.

1. Introduction

Industrial and agricultural activities have been disturbing the global biogeochemical cycles of N, P and C, threatening both ecosystems and the present and future safety and feasibility of humanity [1,2,3]. The increasing availability of P and biologically reactive N in aquatic ecosystems has led to their ecological deterioration globally [4,5] and has compromised their water quality [6,7]. Anoxic conditions are among the most critical effects of cultural eutrophication because they can affect the budgets of N, P and C, and in particular, they can further significantly modify the N transformations and dynamics [8].

Nevertheless, lakes and reservoirs also have the potential to act as N, P and C removers [9], by burying nutrients and organic matter into the sediments [5,10] or by permanently eliminating it through the production of gaseous species such as N₂ and N₂O [9,11] via denitrification [12], anammox and denitrifying anaerobic methane oxidation (DAMO) [13,14]. It has been estimated that 19.7 × 10⁹ kg N y⁻¹ are removed every year globally from watersheds by lentic systems [15]; an important fraction of which (33%) is thought to take place in small reservoirs (surface area between 0.001–50 km²; [15,16] despite that they represent only 6% of the global lentic surface area [17]. This is attributed to their: (1) long residence times compared to rivers and streams, (2) high ratio of surface catchment area−lake surface area, (3) high apparent setting velocities for N and (4) elevated average N loads [15,18]. This removal is enhanced in systems with conditions favorable for denitrification, such as hypoxic or anoxic bottom water and/or sediments, and with a high content of organic carbon and nitrate [18,19].

The importance of reservoirs in N cycling has been increasing globally over the last 60 years [20], not only at the local scale, but also at the regional and global scales [15]. Furthermore, reservoirs and lakes exhibit important spatial variations of the N removal potential [19] and the regulation mechanisms need to be further explored in detail. The need to document N removal variations and the mechanisms that affect them is paramount in tropical latitudes where aquatic systems are highly pressured by anthropic activities [21,22,23]. Because of the higher loads they receive, understanding the N dynamics of tropical reservoirs has become a priority.

Additionally, many reservoirs are also important water sources for densely populated areas and as a result, important water level fluctuations (WLF) take place in these systems. Water depth has been identified as a key factor that underpins the relative roles and fates of N and P in lakes [24]. It has also been suggested that WLF might significantly affect the conditions that regulate N cycling and, therefore, also its removal [18,19]. Such is the case of Valle de Bravo (VB) lake, which supplies Mexico City and other metropolitan areas in the Mexican highlands. Because of this, the basic limnological aspects of VB are reasonably well known [25,26,27,28,29,30,31,32], so it is a convenient system for the exploration of N dynamics. In an initial (2002–2005) time-series mass balance assessment in VB during which high WLF did not occur, Ramírez-Zierold et al. [29] estimated that N₂ fixation exceeded denitrification and doubled the N loading that the reservoir received from rivers and sewage at the time. Nevertheless, important WLF have occurred since then in VB and their significant effect on vertical mixing has been outlined [33].

In the present study we extend and improve the N and P mass balances along a long-term series (16 years) from 2003 to 2018, during which important WLF of up to 12 m took place in the VB reservoir. Our goal is to assess N and P retention, N and P net sedimentation and N net atmospheric fluxes at an ecosystem scale and document their variability as a function of the forcing factors involved, including external N and P loading, limnological conditions and water level fluctuations. We expect that our analysis of this long-term monitoring and mass balances will shed light on the relative relevance of external load variations and water level fluctuations on these fluxes and on the main N cycling processes to contribute to extending the information on tropical reservoirs and to increase the insight on the factors affecting N and P cycling in aquatic ecosystems.

2. Materials and Methods

2.1. Study Area: Valle de Bravo

Valle de Bravo (VB) lake (19°21′30′′ N, 100°11′00′′ W) is the biggest reservoir of the Cutzamala System (CS), formed by seven reservoirs located in different basins (Figure 1), which supplies water to Mexico City and its metropolitan area to the benefit of 4.1 million inhabitants [34]. The VB reservoir has a maximum capacity of 395.5 × 10³ m³, a mean surface area of 18.55 × 10⁶ m² and a mean depth of 21.1 m. However, when intense level fluctuations (1817.5 to 1830.0 m a.s.l) take place, its surface area can decrease down to 13.2 × 10⁶ m². When it is required, water can be pumped back to VB from the reservoirs on the western part of the Cutzamala System (see [29] for details).

Like many other reservoirs in Mexico, Valle de Bravo (VB) shows a detriment in water quality mainly due to its eutrophication since 1992 [35] and health risks due to the presence of cyanotoxins like microcystin and anatoxin-a [31,36,37,38]. For these reasons, information on nutrients (particularly N and P cycling) in this system has become a priority in recent years.

Water arrives to VB through five tributary rivers and three domestic sewage discharges (Figure 2), among which the Amanalco river is the main source [29]. Agriculture in the Amanalco watershed drives increased N and P loads derived from the overuse of fertilizers. Upstream from the Amanalco river, 94 trout aquaculture facilities operate and only 30% of them dispose of their wastewater under Mexican legislation [39]. All of the above have been identified as causes of the eutrophication of VB [26,29].

VB is a warm monomythic water body with strong WLF linked to water management for supply purposes [33]. VB is located in a mountainous region where strong diurnal wind flows [26]. Both WLF and strong wind flows impact diverse limnological and ecological aspects such as stability of stratification, nutrient exchange between the epilimnion and the hypolimnion [33], composition of the phytoplankton and zooplankton community and ecosystem metabolism [31,32,40].

In their preliminary mass balance assessment for the 2002–2005 period, Ramírez-Zierold et al. [29] found that N₂ fixation exceeded denitrification of VB and that the reservoir’s net N₂ fixation doubled the N loading it received from the rivers and sewage together. However, in the following years, important WLF occurred in VB. Fortunately, the continuous monitoring of the reservoir maintained since then makes possible the assessment of the effects of WLF on N and P fluxes, as was previously done for vertical water fluxes [33].

2.2. Samples and Data

The database was comprised of monthly samplings from January 2003 to December 2018, during which the VB reservoir, five rivers (Amanalco, Molino, González, Carrizal and Santa Monica) and three sewage discharges (Tizates, embarcadero I and embarcadero II) were monitored every 30 ± 1 days (N = 200). Flow was measured at rivers and discharges using a drifting buoy and cross-section determinations (cf. [29] for details). For the reservoir, a single central station was monitored after 2006, as [31,32,41] found it to be representative of the reservoir because of its horizontal homogeneity.

Water samples were collected with a 1.5 L Uwitec sampling bottle at all rivers and sewage discharges as well as at the central station of VB at depths of 0, 1, 2, 4, 8, 12, 16, 20, 24 m and the bottom. Samples for total nitrogen (TN) and total phosphorous (TP) determinations were collected in duplicate in 30 mL polypropylene bottles. Using GF/F Whatman^® filters, 60 mL of water were filtered and the filters were stored for particulate organic nitrogen (PON) and particulate organic phosphorous (POP) determinations. All samples were maintained at 4 °C in the dark until the analysis. Additionally, at the central station, vertical profiles of temperature and dissolved oxygen (DO) were done using a Hydrolab DS4/SVR4 field probe from 2003 to 2009 and a multiparametric sonde YSI 6600 from 2009 to 2018.

2.3. Nutrient Analysis

Water samples for the determination of TN and TP, as well as filters for PON and POP, were subjected to high temperature persulfate oxidation to oxidize all the P and N species to nitrate and soluble reactive phosphorous (SRP), following [42]. Nutrient (SRP, ammonia, nitrite and nitrate) concentrations were determined using a Skalar San-Plus segmented-flow analyzer [43,44]. Dissolved inorganic nitrogen (DIN) was calculated as ammonia + nitrite + nitrate.

2.4. Mass Balance Approach

The water mass balance (Equation (1)) considered the reservoir as a single box and assumed that groundwater seepage was negligible [31].

$$\mathit{\Delta }V/\mathit{\Delta }t=\Sigma R_i+P_i-W-\Sigma \left(E-P\right)\times A$$

(1)

where $\raisebox{1ex}{$\mathit{\Delta }V$}\!\left/ \!\raisebox{-1ex}{$\mathit{\Delta }t$}\right.$ is the volume change of the reservoir, R_i is the sum of all the river and sewage inflow rates, P_i is the pump-back inflow, W is the water withdrawal rate, E is evaporation, P is precipitation (m³d⁻¹ in all cases) and A is the lake surface area (m²).

While the volume change of the reservoir was available on a daily basis (provided by CONAGUA, National Water Commission), the rivers and sewage inflows could only be measured once a month since no government monitoring of the rivers was done at this time. Therefore, to improve the representativeness of our river discharge assessment, we calculated the water residuals (Wres, Equation (2)), found using our monthly measurements and assuming a closed water mass balance (with no other significant inputs or outputs). These residuals were then used to correct the inflows $\left({\overline{R}}_i\right)$ by iteration (3), redistributing the $Wres$ water flow among the rivers according to each river’s contribution fraction (x_i, ranging between 0 and 1) to the total river input, which were assessed using the river´s average contribution fraction during the years for which the minimal deviations $\left(Wres\right)$ were observed for the annual mass balances.

$$Wres=\raisebox{1ex}{$\mathit{\Delta }V$}\!\left/ \!\raisebox{-1ex}{$\mathit{\Delta }t$}\right.-\left(\Sigma R_i+P_i-W-\Sigma \left(E-P\right)\times A\right)$$

(2)

$${\overline{R}}_i=R_i+x_i\left(Wres\right)$$

(3)

$$\raisebox{1ex}{$\mathit{\Delta }V$}\!\left/ \!\raisebox{-1ex}{$\mathit{\Delta }t$}\right.=\Sigma {\overline{R}}_i+P_i-W-\Sigma \left(E-P\right)\times A$$

(4)

2.5. Estimation of N and P Net Internal Fluxes

The phosphorous and nitrogen mass balances were represented by Equation (5), as proposed by Ramírez-Zierold et al. [29]:

$$\raisebox{1ex}{$\mathit{\Delta }{M}_{\left(P,N\right)}$}\!\left/ \!\raisebox{-1ex}{$\mathit{\Delta }t$}\right.=\Sigma I_{\left(P,N\right)}-O_{\left(P,N\right)}\pm NI{P}_{\left(P,N\right)}$$

(5)

where $\raisebox{1ex}{$\mathit{\Delta }{M}_{\left(P,N\right)}$}\!\left/ \!\raisebox{-1ex}{$\mathit{\Delta }t$}\right.$ is the mass change of N or P in the reservoir between consecutive sampling surveys. P and N inputs ($I_{\left(P,N\right)})$ were calculated using the mean concentrations of TP and TN determined on each source during both samplings multiplied by the corrected inflow $\left({\overline{R}}_i\right)$. Output ($O_{\left(P,N\right)})$, was calculated using the TP and TN mean concentration near the bottom of the reservoir lake (from where the water is withdrawn, [29]) multiplied by the water withdrawal rate.

$NI{P}_{\left(P,N\right)}$ represents the net element flux due to all the biogeochemical processes that take place internally within the reservoir. Because there are no processes that drive an exchange with the atmosphere, in the case of P, $NI{P}_{\left(P\right)}$ is the net exchange between the water column and the sediments. For N, nevertheless, $NI{P}_{\left(N\right)}$ includes the net N exchange with the sediments and also the exchange between the water column and the atmosphere (Equation (6)) due to the input of N₂ fixation and the N emissions resulting from multiple processes such as denitrification, anammox and DAMO.

To estimate the N net flux between the water column and the sediments $\left(NI{P}_{\left(N\right)Sed}\right)$, we used the P net sedimentation ($NI{P}_{\left(P\right)})$ and multiplied it by the particle N:P ratio (PON:POP) found in the hypolimnion of VB, as this is the lake’s layer from which both sedimentation and exchange directly take place. This meant a refinement of the $\left(NI{P}_{\left(N\right)Sed}\right)$ estimation, as compared to that of [29], where the full water column N:P ratio was used.

The net N flux between the atmosphere and the water column was the assessed from the $NI{P}_{\left(N\right)}$ after subtracting $NI{P}_{\left(N\right)Sed}$ from it in Equation (6):

$$NI{P}_{\left(N\right)}=NI{P}_{\left(N\right)Atm}+NI{P}_{\left(N\right)Sed}$$

(6)

2.6. Data Statistical Analysis and Data Visualization

Pearson’s correlations were used to find significant relationships between net N fluxes and water level, environmental parameters and vertical fluxes of DIN [33]. Data rearrangement and their subsequent analyses were done using Microsoft Excel (King County, WA, USA), Ggplot2 [45] and Tidyverse [46] packages for R software [47]. To visualize and assess interannual and long-term variations, the data were integrated on a yearly basis. To help visualize the seasonal variations, the $NI{P}_{\left(N\right)Atm}$ data series was smoothed with a three-phase moving average until a clear seasonal pattern variation was revealed. Then, the data were averaged on a monthly scale to obtain a mean annual cycle of $NI{P}_{\left(N\right)Atm}$ on which the relative importance of the different N processes could be assessed.

3. Results

3.1. Water Budget

Throughout the whole 2003–2018 period, the main water inflow to VB came from the rivers, mainly Amanalco and Molino, and the outflow was due to withdrawal (Figure 3). During the dry seasons (from November to May), the rivers decreased their inflows while withdrawal increased to provide water to Mexico and other cities. During the rainy season, the rivers’ input increased, and the water withdrawal decreased. This caused a discrete seasonal oscillation of the reservoir’s water level, in the range of 4–6 m. However, as discussed by [35], during some years, the water imbalance was higher, causing higher WLF of up to 12 m below the reservoir’s capacity.

Two statistically different (t(9) = 4.75, p = 0.001) water level patterns were identified for 2003–2018: (1) low water level fluctuation (L-WLF) periods when water level remained within 4 to 6 m of the capacity level (1830.0 m a.s.l), which occurred from 2003 to 2005 (L-WLF-1) and from 2014 to 2018 (L-WLF-2) (M = 1827.4, SD = 0.81) and (2) a high water level fluctuation (H-WLF) period that occurred from 2006 to 2013 when the reservoir lake’s water level decreased farther below the historic mean (M = 1824.0, SD = 1.88). Water budget means for each of the periods can be found in Table S1 (Supplementary Material). During the H-WLF period, the input of rivers and sewages decreased ~16% overall compared to the L-WLF periods and the water level reached extraordinary minima (<1820 m a.s.l) in 2006, 2009 and 2013 (Figure 3). Water management during this period implied a decrease (by 21%) of withdrawal and an increase in pump-back from other reservoirs (from <4% to >10% of the inputs) in order to recover the water level. As a result, a recovery of the reservoir water level to its average level was achieved in ~5 months, particularly between 2009 and 2013.

3.2. N and P Mass Budgets

As previously found by Ramírez-Zierold et al. [29], both the N and P budgets and the resulting NIPs varied significantly at the monthly scale. Therefore, to smooth out the short-term variability and to facilitate the visualization of the long-term trends as well as the net differences among years and among different WLF periods (

Table 1. Phosphorous and nitrogen budgets for Valle de Bravo during 2003–2018. Mean annual phosphorous (TP) and nitrogen (TN) budget components of Valle de Bravo are shown for: the low water level fluctuation-1 (L-WLF-1) period (2003–2005), the high water level fluctuation (H-WLF) period (2006–2013) (shaded in purple for better visualization), the water level fluctuation-2 (L-WLF-2) period (2014–2018), the overall means for the full study period (2003–2018) as well as the percentage of each flux relative to the total input. The net internal fluxes of P and N are shaded in orange, the net sedimentation fluxes are shaded in yellow, and the atmosphere−water fluxes of N are shaded in blue.

	L-WLF-1	H-WLF	L-WLF-2	Overall
TP	(2003–2005)	(2006–2013)	(2014–2018)	(2003–2018)
Mean mass reservoir (103 kg P)	26.3	30.8	26.5	28.6
Δ mass (103 kg P y−1)	1.3	2.7	1.7	2.1
Fluxes (103 kg P y−1)					% Input
Rivers	60.7	56.4	34.0	50.2	73%
Sewage	23.2	15.9	7.5	14.7	21%
Pump-back	1.1	4.7	3.7	3.7	5%
Total input	85.0	77.1	45.2	68.6	100%
Output	−21.0	−20.5	−17.4	−19.6	−29%
Input−Output	64.0	56.6	27.8	49.0	71%
Net internal flux	−62.7	−53.9	−26.1	−46.9	−68%
Net sedimentation	−62.7	−53.9	−26.1	−46.9	−68%
TN					% Input
Mean mass reservoir (103 kg N)	382.6	289.0	288.6	306.4
Δ mass (103 kg N y−1)	65.1	−22.7	−0.3	0.8
Fluxes (103 kg N y−1)
Rivers	354.6	293.3	281.0	300.9	80%
Sewage	51.7	58.5	46.0	53.3	14%
Pump-back	8.0	27.7	26.0	23.5	6%
Total input	414.3	379.5	353.0	377.7	100%
Output	−286.4	−194.0	−194.7	−211.6	−56%
Input−Output	127.9	185.4	158.3	166.2	44%
Net internal flux	−62.8	−208.2	−158.5	−165.4	−44%
NOP: POPhipo	36	26	33	29.9	8%
Net sedimentation	−1049.9	−662.4	−497.4	−683.5	−181%
Atmosphere−water fluxes	987.2	454.2	338.8	518.1	137%

), we integrated the N and P budgets on the annual scale (Figure 4 and Figure 5). The H-WLF period is shaded in purple (as in Table 1 and Table S1, Supplementary Material) to aid in the identification of the high WLF.

P budget (Figure 4) shows that VB received an average P load of 68.6 × 10³ kg P y⁻¹ during 2003–2018. Most of the P input to VB (73% of the overall mean) arrived through the rivers but the input from sewages and pump-back from other reservoirs of the Cutzamala System was also significant. The withdrawal of water from the reservoir removed an average of 19.6 × 10³ kg P y⁻¹ during the same period, only about one third (29% of the overall mean) of the P arriving to VB, so the rest of the P load (68% of the overall mean), would have been incorporated into the sediments as a result of the net P sedimentation flux $NI{P}_{\left(P\right)}$ (Figure 4), which averaged 46.9 × 10³ kg P y⁻¹ for the whole sampling period.

Overall, P inputs to VB decreased throughout 2003–2018 at an average rate of −3.4 × 10³ kg P y⁻¹. As a result, and because the amount removed from the lake by water withdrawal remained nearly constant, the amount of P removed into the sediments also decreased throughout the period at a mean rate of −2.8 × 10³ kg P y⁻¹. The yearly variations observed in the inputs were reflected in the net P sedimentation, with 2010 having the maximum P input and sedimentation annual rates and 2018 having the lowest rates, when the P budget even yielded a net P efflux of P from the sediments to the water column on an annual basis. The H-WLF period did not show a pronounced difference when compared to the L-WLF periods and the mean annual fluxes were within 10% of the mean values for the previous L-WLF-1 (Table 1). Nevertheless, after the H-WLF period, the P net sedimentation flux was on average 26.1 × 10³ kg P y⁻¹, less than half of those occurring during the two previous periods, either under low WLF (62.7 × 10³ kg P y⁻¹, L-WLF-1) or high WLF (53.9 × 10³ kg P y⁻¹, H-WLF).

N loading to VB was also high during 2003–2018 when it received an average of 377.7 × 10³ kg N y⁻¹ through the water inputs to the reservoir (Table 1). The main nutrient input was from rivers (80% overall) and sewages (14% overall). Nevertheless, the increase in the pump-back operations from other reservoirs during the H-WLF period meant almost a four-fold increase in the N input to VB through this source, which passed from being <2% of the total N input during the L-WLF-1 to <7% during the H-WLF period (Figure 5) and remained near this proportion during the rest of the study period (L-WLF-2). In spite of this, N inputs from the water sources showed an decreasing trend (−6.5 × 10³ kg N y⁻¹) overall throughout the study period, likely as a result of the management actions taken on the rivers and sewage inputs [29,38].

N internal ($NI{P}_{\left(N\right)}$) varied significantly during 2003–2018 as compared to P, although VB also retained N (i.e., input > output) throughout 2003–2018, behaving as an N sink (Table 1; Figure 5). VB had a net N flux to the sediments ($NI{P}_{\left(N\right)Sed}$) and a net N flux from the atmosphere ($NI{P}_{\left(N\right)Atm}$) throughout the sampled period but the magnitude of these fluxes changed markedly among the years. Both net fluxes decreased significantly during the H-WLF period, particularly in the lowest level years (2006, 2009 and 2013). As a result, $NI{P}_{\left(N\right)Sed}$ decreased on average by 37% during the H-WLF period relative to average of the previous years (L-WLF-1), and $NI{P}_{\left(N\right)Atm}$ decreased to around almost half (54%) compared to the same years (Table 1). In spite of the return to low level fluctuations after 2013, both net N fluxes remained just as low on average during L-WLF-2 and, in fact, continued decreasing. Overall, $NI{P}_{\left(N\right)Sed}$ decreased on average by 35 × 10³ kg N y⁻¹ and $NI{P}_{\left(N\right)Atm}$ by 30 × 10³ kg N y⁻¹ throughout the 2003–2018 period, an impressive decrease of about 3% per year in the nitrogen fluxes of the system (Figure 5).

Correlation analysis (Figure 6) shows that the variations of the net N atmospheric flux $NI{P}_{\left(N\right)Atm}$, were closely related to the water level changes that occurred in VB during the H-WLF period. $NI{P}_{\left(N\right)Atm}$ decreased significantly as the water level drooped and the stability of the water column diminished. It also correlated significantly, but inversely, with the oxygen concentration and the DIN vertical flux, which increased as the level dropped and vertical mixing between the epilimnion and the hypolimnion was intensified [35], particularly during the H-WLF years.

Within the annual scale, after the smoothing of short-term variations, the net N atmospheric flux $NI{P}_{\left(N\right)Atm}$ showed a clear seasonal pattern associated with the annual stratification−circulation cycle of this monomictic lake (Figure 7). Overall, $NI{P}_{\left(N\right)Atm}$ increased during the stratification period when the conditions are likely most favorable for N₂ fixation, peaking during June. In contrast, $NI{P}_{\left(N\right)Atm}$ was lowest during the circulation period, particularly during its second half in the first months of this year. At this time, nitrate peaks in VB, so when stratification begins, denitrification is likely favored, delaying the evident increase of $NI{P}_{\left(N\right)Atm}$ until May when the stratification is clearly established and nitrate has been depleted. As the cycle proceeds, small observed decrements in $NI{P}_{\left(N\right)Atm}$ (as observed in April, July and September) may be due to other processes that might cause the efflux of N from the system, such as Annamox and DAMO, as suggested by unpublished data [48]. Finally, because ammonia is accumulated in the epilimnion during the stratification period and $NI{P}_{\left(N\right)Atm}$ decreases sharply as its annual circulation begins, it is posed that significant emissions of ammonia [49] may occur at this time of the year in VB.

Independently of these variations, the mean values of all the components of the N budget are useful to make a graphical summary and compare their relative magnitudes (Figure 8). As done in Table 1 and Table S1 (Supplementary Materials) and in previous assessments (i.e., [29]), we compared the N fluxes with the external loading to the system. It is interesting to note that on average, VB transferred to its sediments 683.5 × 10³ kg N y⁻¹ (181 % of the external load it received), while 211.6 × 10³ kg N y⁻¹ (56 % of this load) was removed via water withdrawal. The apparent deficit of N could suggest this is caused by the high net N influx from the atmosphere (518.1 × 10³ kg N y⁻¹ on average), which is about 137 % higher than the external N loading (Figure 8).

4. Discussion

Sustained monitoring of VB lake showed that during the 16-year period of 2003–2018, this tropical reservoir of the CS continued receiving important loads of P and N from human activities in its basin, mainly through rivers and sewages, but also, occasionally, through pump-back operations from other reservoirs of the CS. These loads, which averaged 68.6 × 10³ kg P y⁻¹ and 377.7 × 10³ kg N y⁻¹, represented 239% in the case of P and 123% for N of the total mass of these nutrients contained in the water of the lake (28.6 × 10³ kg P y⁻¹ and 306.4 × 10³ kg N y⁻¹, respectively) on average during this period. So, on a yearly basis, the system received almost 2.5 times the amount of P it already contained and 1.2 times the N that was already there. These proportions help visualize the magnitude of the environmental pressure VB has to handle as an ecosystem and as a biogeochemical processor.

4.1. Most of the External Load Translates into Net Sedimentation

The main consequence of the external load of N and P to VB—besides eutrophication and all its associated impacts and risks [26]—is a high net sedimentation flux of P and N to the sediments of the lake, which we estimate is the destiny of 68% of the P load and is equivalent to 44% of the N load on average. The rest of the external load to the reservoir is removed by water withdrawal and so 29% of the P load and 56% of the N load are transferred to the water supply processing system of the CS.

Although these loads are still very high, they likely evidence an improvement as compared to the loads of 120.8 × 10³ kg P y⁻¹ and 591.8 × 10³ kg N y⁻¹ assessed by Ramírez-Zierold et al. [29] in a preliminary assessment of 2002–2005 conditions. The detailed monitoring of P and N that has proceeded since then allowed the verification of a decreasing trend of external inputs of P (~5% per year) and N (~2%) to VB for the 2003–2018 period overall. This trend was particularly evident after 2010 when the preliminary assessment was published and its results promoted the actions of emerging non-governmental organizations like Pro-Valle (e.g., [50] and the onset of citizen science actions in the basin, which in turn spurred a series of management actions directed to reduce the external loads to VB), included waste water treatment from agriculture, livestock and aquaculture activities on the Amanalco and Tizates river basins [38,51,52].

Our results also show that a sustained decrease of the net sedimentation fluxes occurred in a concomitant way and a trend parallel to the decrease of the loads, particularly in the case of P. In this case, the average decrease rate is ~6% of the net P sedimentation rate, relatively slightly higher than the average P load decrease rate. This is an encouraging finding that is backed up by the consistency of our estimates of net P sedimentation with external P loading. Our mean net P sedimentation of 2.7 g P m⁻² y⁻¹ is quite close to the average of 3.3 g P m⁻² y⁻¹ reported for tropical reservoirs [22] and for other aquatic systems such as Lake Kariba [53] and Lake Xuanwu [54]. At the same time, our estimates for the initial 2003–2005 period also match very well with the P sedimentation rates in VB estimated for those years (2.0 to 4.2 g m⁻² y⁻¹) by Carnero-Bravo et al. [55], who used radiochronology studies to assess sedimentation in VB from 1992 to 2006. These results also sustain that the net P sedimentation preliminary assessed by Ramírez-Zierold et al. [29] for 2002–2006 (~6 g m⁻² y⁻¹) was likely an overestimation, as suspected.

Similarly, in the case of the net N sedimentation rate, our estimates for 2003–2005 also indicate an overestimation in the preliminary assessment (75 g m⁻² y⁻¹) of [29] for the period of 2002–2005 and yield a much lower overall mean (41.9 g m⁻² y⁻¹) for the whole 2003–2018 period. This difference is likely partially due to improvements in the assessment approach made here, both on the water balance (i.e., the residual redistribution on the river inputs) and on the use of the hypolimnetic PON:POP proportion instead of the proportion within the entire water column. This second modification may be relevant because important groups of the epilimnetic phytoplankton biomass, like cyanobacteria, might not reach the sediments as much as others, like diatoms. These new estimates are much closer to the N sedimentation estimated by radiochronology [55], which ranged from 4.9 g m⁻² y⁻¹ as a basal value before 1992 to 48 g m⁻² y⁻¹ in 2005, solving much of the suspected overestimation problem. The direct measurements of sedimentation rates could help to further confirm these results.

4.2. Net N Atmospheric Flux Was Lower during High Water Level Fluctuations

As previously described, net N atmospheric flux decreased drastically (by 54%) and significantly during the high WLF period as compared to the previous low WLF period (Table S1, Supplementary Materials). Furthermore, within the H-WLF, net N atmospheric flux correlated significantly with the water level change over the years (n = 9, r² = 0.58; Figure 6a). In a detailed analysis of the effects of high WLF in VB, Merino-Ibarra et al. [33] found that during the high WLF years, as the water level decreased, vertical mixing between the epilimnion and the hypolimnion of VB intensified through boundary mixing events without breaking the stratification. This enhanced mixing would imply an important vertical exchange of nutrients and oxygen between these two layers.

The input of oxygen to the otherwise anoxic hypolimnion and the vertical input of nutrients to the impoverished epilimnion could certainly affect the multiple processes involved in N cycling in diverse ways. We can identify three mechanisms that may be occurring in VB when vertical mixing intensifies and the net N atmospheric flux decreases.

One mechanism would be the direct inhibition of N₂ fixation due to the additional input of DIN to the epilimnion (Figure 6b). When the availability of inorganic nitrogen species in the epilimnion is high, N₂ fixation decreases, as cyanobacteria (e.g., Nostocales, found in VB [31]) suppress the synthesis of the new nitrogenase to reduce the high energetic cost involved in its production and therefore, their N₂ fixation potential is diminished. Another mechanism driving a decrease in N₂ fixation could be the effect of mixing on the abundance of N fixers. It has been demonstrated that an increase in turbulence in the water column plays a profound impact on phytoplankton composition, including the promotion of non-cyanobacterial phytoplankton (e.g., diatoms and chlorophyta) [31]. In VB, diatoms such as Cyclotella ocellata, Fragilaria crotonensis were abundant during the low-level stratifications of 2008 and 2009, while Cyanophyta was notoriously abundant during high-level 2001–2002, 2010 and 2017 stratifications [28,36,37,56,57]. Additionally, zooplankton community composition changes (i.e., increasing abundance of cladoceran and copepod assemblages during low-level stratifications instead of rotifer-dominated zooplankton communities) could also impact phytoplankton through herbivory but also sedimentation, primary and secondary production [30].

The other way mixing could have decreased the net N atmospheric flux is through favoring or enhancing the diverse processes that can cause N loss to the atmosphere. The exchange of oxygen, ammonium and nitrate between the oxic and anoxic layers during stratification creates suitable conditions for nitrification coupled with denitrification and other microbial processes, promoting N loss from the water column to atmosphere, as observed in several aquatic ecosystems [9,14,58,59]. Denitrification strongly depends on the nitrate availability in the water column [9,19,60], which could increase in association with the inputs of oxygenated water to the hypolimnion of VB, which has a high ammonium content, favoring denitrification episodes. In this study, using a mass balance approach, we cannot estimate the denitrification rate itself, but the inverse correlation between the net N atmospheric flux and the oxygen content in the water column (Figure 6d), as well as its direct correlation with the mean nitrate concentration in the water column (not shown) supports this possibility in VB.

Although the net sedimentation fluxes also diminished from the first low WLF period (L-WLF-I) to the high WLF period (H-WLF), it was much lower than that of net N atmospheric flux (54%), particularly for net P sedimentation (only 14%). This decoupling of the two fluxes could be explained by the enriched abundance of diatoms as a response to the turbulence rising during H-WLF [31]. Diatoms have been cited to exhibit higher sedimentation rates than other phytoplankton species [61,62]. In contrast, diverse species of cyanobacteria [63], in particular some of the genus Microcystis, are known to possess gas vacuoles to regulate their buoyancy [64]. Additionally, cyanobacterial taxa are known to produce mucilage or polysaccharide-rich extracellular polymeric substances (EPS) that also contribute to buoyancy [65,66]. So, the decrease of cyanobacteria dominance during the H-WLF period would diminish the part of the biomass that is mineralized in the water column rather than settling in the lake sediment [63] and would therefore also contribute to a relatively smaller decrease of net sedimentation. The higher decrease in net N sedimentation (47%) as compared to the P decrease (14%) could be due in part to the relative impoverishment of N in the planktonic community concomitant with the decreased N₂ fixation during this period.

Overall, our results also raise interesting possibilities and questions that might deserve further analysis and complementary investigations. For example, both the net sedimentation and net atmospheric fluxes were relatively low during the second low WLF period (L-WLF-2) as compared to the first low WLF, suggesting that some of the changes that occurred during the H-WLF period may have prevailed or may have modified the ecosystem in a lasting way. Alternatively, this may be the result of the long-term decrease in the external loads derived from the social demand of management actions.

Our results in VB suggest that strong water level fluctuations affect, not only their temperature, oxygen, nutrients, and planktonic composition, but they also impinge on N cycling, particularly by decreasing the net N atmospheric flux in the water column, so they are a potentially valuable management tool, and should be considered in water management plans to relieve eutrophic conditions in the reservoir. Further studies addressed to assess specifically N₂ fixation, nitrification, denitrification and other N processes will be useful to further examine our preliminary conclusions.

5. Conclusions

During the 16 years of this long-term mass balance study, the tropical reservoir lake of Valle de Bravo was a sink of N (−41.7 g N m⁻² y⁻¹) and of P (−2.7 g P m⁻² y⁻¹), mainly due to a net sedimentation flux.

The P retention through net sedimentation flux in this system was driven mainly by the external load to the reservoir, while in the case of N, a very high net atmospheric flux of 31.6 g N m⁻² y⁻¹, on average, was an additional driver contributing to the high retention of this element.

Overall, the high N net atmospheric exchange found in VB was likely dominated by the N input due to N₂ fixation over the combined output of the multiple processes that drive N loss, such as denitrification, anammox and DAMO.

The long-term database used allowed assessing the effects of water level fluctuations in tropical reservoirs, revealing that during high water level fluctuation periods, the N atmospheric flux decreased by half.

We propose that during low water level years, vertical mixing boundary processes likely enhance denitrification through the oxygenation of the water column and inhibit N₂ fixation through the injection into the epilimnion of inorganic N from the hypolimnion.

These results show that water level regulation can be used as a useful management tool to enhance N retention in aquatic ecosystems, through decreasing N₂ fixation and promoting denitrification and other microbial processes that can help alleviate the N load to these systems.