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

: 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 − 2 y − 1 ) and P ( − 2.7 g P m − 2 y − 1 ), mainly occurring through net sedimentation, equivalent to 181% and 68% of their respective loading (23.0 g N m − 2 y − 1 and 4.2 g P m − 2 y − 1 ). The N mass balance also showed that the Valle de Bravo reservoir has a high net N atmospheric inﬂux (31.6 g N m − 2 y − 1 ), which was 1.3 times the external load and likely dominated by N 2 ﬁxation. P ﬂux was driven mainly by external load, while in the case of N, net ﬁxation also contributed. During a period of high water level ﬂuctuations, the net N atmospheric ﬂux 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

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)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(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 −2 y −1 ) and P (−2.7 g P m −2 y −1 ), mainly occurring through net sedimentation, equivalent to 181% and 68% of their respective loading (23.0 g N m −2 y −1 and 4.2 g P m −2 y −1 ). The N mass balance also showed that the Valle de Bravo reservoir has a high net N atmospheric influx (31.6 g N m −2 y −1 ), which was 1.3 times the external load and likely dominated by N 2 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.

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 2 and N 2 O [9,11] via denitrification [12], anammox and denitrifying anaerobic methane oxidation (DAMO) [13,14]. It has been estimated that 19.7 × 10 9 kg N y −1 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 2 ; [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)(2003)(2004)(2005) time-series mass balance assessment in VB during which high WLF did not occur, Ramírez-Zierold et al. [29] estimated that N 2 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 longterm 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.

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 3 m 3 , a mean surface area of 18.55 × 10 6 m 2 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 6 m 2 . 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). 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.

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 3 m 3 , a mean surface area of 18.55 × 10 6 m 2 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 6 m 2 . 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]. 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 N2 fixation exceeded denitrification of VB and that the reservoir's net N2 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 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 2 fixation exceeded denitrification of VB and that the reservoir's net N 2 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].

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.

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.

Mass Balance Approach
The water mass balance (Equation (1)) considered the reservoir as a single box and assumed that groundwater seepage was negligible [31].
where ∆V/∆t 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 3 d −1 in all cases) and A is the lake surface area (m 2 ). 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 R i 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 (Wres) were observed for the annual mass balances.

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]: where ∆M (P,N) /∆t is the mass change of N or P in the reservoir between consecutive sampling surveys. P and N inputs (I (P,N) ) were calculated using the mean concentrations of TP and TN determined on each source during both samplings multiplied by the corrected inflow R i . Output (O (P,N) ), 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. N IP (P,N) 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, N IP (P) is the net exchange between the water column and the sediments. For N, nevertheless, N IP (N) 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 2 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 N IP (N)Sed , we used the P net sedimentation (N IP (P) ) 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 N IP (N)Sed 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 N IP (N) after subtracting N IP (N)Sed from it in Equation (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 N IP (N)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 N IP (N)Atm on which the relative importance of the different N processes could be assessed.

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) 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.
Water 2022, 14, x FOR PEER REVIEW 7 of 19 discussed by [35], during some years, the water imbalance was higher, causing higher WLF of up to 12 m below the reservoir´s capacity.  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.

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), we integrated the N and P budgets on the annual scale (Figures 4 and 5). The H-WLF period is shaded in purple (as in Tables 1 and 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 3 kg P y −1 during 2003-2018. Most of the P input to VB (73% of the overall mean) arrived through

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), we integrated the N and P budgets on the annual scale (Figures 4 and 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 3 kg P y −1 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 3 kg P y −1 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 N IP (P) (Figure 4), which averaged 46.9 × 10 3 kg P y −1 for the whole sampling period.   Overall, P inputs to VB decreased throughout 2003-2018 at an average rate of −3.4 × 10 3 kg P y −1 . 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 3 kg P y −1 . 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 3 kg P y −1 , less than half of those occurring during the two previous periods, either under low WLF (62.7 × 10 3 kg P y −1 , L-WLF-1) or high WLF (53.9 × 10 3 kg P y −1 , H-WLF).
N loading to VB was also high during 2003-2018 when it received an average of 377.7 × 10 3 kg N y −1 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 3 kg N y −1 ) overall throughout the study period, likely as a result of the management actions taken on the rivers and sewage inputs [29,38].
N internal (N IP (N) ) 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 (N IP (N)Sed ) and a net N flux from the atmosphere (N IP (N)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, N IP (N)Sed decreased on average by 37% during the H-WLF period relative to average of the previous years (L-WLF-1), and N IP (N)Atm decreased to around almost half (54%) compared to the same years (Table 1  Overall, P inputs to VB decreased throughout 2003-2018 at an average rate of −3.4 × 10 3 kg P y −1 . 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 3 kg P y −1 . 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  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 3 kg N y −1 ) overall throughout the study period, likely as a result of the management actions taken on the rivers and sewage inputs [29,38].
N internal ( ( ) ) 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 ( Correlation analysis ( Figure 6) shows that the variations of the net N atmospheric flux N IP (N)Atm , were closely related to the water level changes that occurred in VB during the H-WLF period. N IP (N)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 N IP (N)Atm showed a clear seasonal pattern associated with the annual stratification−circulation cycle of this monomictic lake (Figure 7). Overall, N IP (N)Atm increased during the stratification period when the conditions are likely most favorable for N 2 fixation, peaking during June. In contrast, N IP (N)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 N IP (N)Atm until May when the stratification is clearly established and nitrate has been depleted. As the cycle proceeds, small observed decrements in N IP (N)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 N IP (N)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.
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 ( ) , were closely related to the water level changes that occurred in VB during the H-WLF period.
( ) 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. increased during the stratification period when the conditions are likely most favorable for N2 fixation, peaking during June. In contrast, ( ) 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 ( ) until May when the stratification is clearly established and nitrate has been depleted. As the cycle proceeds, small observed decrements in ( ) (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 ( ) 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 Tables 1 and 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 3 kg N y −1 (181 % of the external load it received), while 211.6 × 10 3 kg N y −1 (56 % of this load) was 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 Tables 1 and 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 3 kg N y −1 (181 % of the external load it received), while 211.6 × 10 3 kg N y −1 (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 3 kg N y −1 on average), which is about 137 % higher than the external N loading (Figure 8).

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 3 kg P y −1 and 377.7 × 10 3 kg N y −1 , 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 3 kg P y −1 and 306.4 × 10 3 kg N y −1 , 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.

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 3 kg P y −1 and 377.7 × 10 3 kg N y −1 , 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 3 kg P y −1 and 306.4 × 10 3 kg N y −1 , 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.

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 3 kg P y −1 and 591.8 × 10 3 kg N y −1 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 −2 y −1 is quite close to the average of 3.3 g P m −2 y −1 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 −2 y −1 ) 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 −2 y −1 ) 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 −2 y −1 ) of [29] for the period of 2002-2005 and yield a much lower overall mean (41.9 g m −2 y −1 ) 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 −2 y −1 as a basal value before 1992 to 48 g m −2 y −1 in 2005, solving much of the suspected overestimation problem. The direct measurements of sedimentation rates could help to further confirm these results.

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 2 = 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 2 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 2 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 2 fixation potential is diminished. Another mechanism driving a decrease in N 2 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 2 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 2 fixation, nitrification, denitrification and other N processes will be useful to further examine our preliminary conclusions.

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 −2 y −1 ) and of P (−2.7 g P m −2 y −1 ), 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 −2 y −1 , 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 2 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 2 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 2 fixation and promoting denitrification and other microbial processes that can help alleviate the N load to these systems.