Modeling the Responses of Dissolved Oxygen and Nitrate Concentrations due to Land Use and Land Cover Change Scenarios in a Large Subtropical Reservoir

_______________________________________________________________ This manuscript is a preprint uploaded to EarthArxiv. This preprint has been submitted for publication in Environmental Modelling & Software and has not yet been peer-reviewed. We welcome feedback, discussion and comments at any time. Feel free to get in touch with one of the authors. _______________________________________________________________


INTRODUCTION
The deterioration of freshwater quality worldwide has resulted in various factors, mainly due to water level changes and inflow nutrient loadings (Gilboa et al., 2014). It is important to identify and quantify the drivers and pressures related to water pollution, especially in urban water bodies (Tuerlinckx et al., 2019). The main ongoing challenges facing inland waters (Downing, 2014) are eutrophication (Liu et al., 2021), agriculture impacts (Lopes et al., 2020), super exploration of the water resource (Simonovic and Arunkumar, 2016), water withdrawal (Feldbauer et al., 2020) and climate change (O'Reilly et al., 2015;Woolway et al., 2020).
According to Gregorio and Jansen (1998), land cover is the observed (bio) physical cover on the earth's surface and land use is explained by human activities, arrangements and inputs to produce, change or maintain some types of land cover. The literature has been investigated internal and external factors associated with consequences on water quality (Laura et al., 2021;Lopes et al., 2020). The main driver associated with the alterations in water quality that causes direct consequences on ecosystem services (e.g., water provisioning for drinking and water purification) is catchments' land use change (Grizzetti et al., 2016). A recent study has proven that total phosphorus release has equal importance to climate change drivers to lead to water quality deterioration in lakes (Shuvo et al., 2021).
Coupled hydrological-biogeochemical models has been used to perform predictions of land use alterations in surrounding watersheds and consequences on water quality (Fenocchi et al., 2020;Messina et al., 2020). Integrating geography information system (GIS) technologies have been giving insightful qualitative answers for water resource management studies (Liu et al., 2018;Soares and Calijuri, 2021).
The Itupararanga reservoir is a large multipurpose reservoir built in the Southeast of Brazil in 1914. The reservoir is classified by its uses as "Class 2 freshwaters" (Brasil, 2005;Melo et al., 2019), which means the drinking water supply, after passing through a treatment plant; the protection of aquatic communities, recreational, irrigational, and fishing activities. The literature has suggested that the main driver for the high trophic state in the Itupararanga reservoir, especially in the riverine zone is the nutrients released from the tributaries Frascareli et al., 2015) and there is a high relationship between nitrogen concentrations and phytoplankton growth, mainly cyanobacteria abundance (Beghelli et al., 2016).
Our hypothesis was that the release of allochthonous inorganic nutrients from the reservoir's tributaries due to land use changes is one of the main drivers of reservoir water quality deterioration. In this study, we aimed to assess the responses of dissolved oxygen and nitrate concentrations of the Itupararanga reservoir applying a biogeochemical model and a simple distributed basin load model. To do this, we evaluated three land use scenarios based on the basin transition potential to the 2050s.

Study site
Itupararanga reservoir is a large lake built in 1914 in the Southeast of Brazil in the Alto Sorocaba basin to support multiple uses, mainly hydropower generation and drinking water supply for almost 1 million people. On the other hand, agriculture is in second place with 40% occupation in the basin. Pasture (12%) and urban areas (3%) followed them.
The main human activities that have compromised the reservoir's water quality are the construction of subdivisions, such as farms and summer houses, intensive use of irrigation and pesticides, and the lack of land use zoning that disciplines the form of disorderly occupation (MANFREDINI, 2018). Another source of pollution in the Itupararanga reservoir is sewage discharge due to poor treatment mainly in Ibiúna, a small city located near the reservoir headwater that releases half of the sewage effluents without any treatment (FABH-SMT, 2018) in the tributary streams.

Database processing
All the processing of the physical variables required by the hydrodynamic model was assumed exactly the same as the previous application of GLM published by Barbosa et al. (2021). Details about the water quality concentration time-series, as well as all the assumptions made, and their processing are given in this section.
Using data from previous studies (Cunha, 2012;Rôdas, 2013;Garcia 2013), it was possible to calculate that the dissolved oxygen (DO) loads in the Sorocaba River upstream of the reservoir were statistically coincident (linear adjustment, r2 =0.85 (n=20)) with the sum of the flows and DO loads of the Sorocabuçu and Sorocamirim rivers measured ~ 3 km from the head of the Sorocaba River. Although upstream Sorocaba River receives effluents from the ETE Ibiúna, this has not sufficiently caused water quality deterioration due to the low flow.
We used flow data along the Sorocabuçu and Sorocamirim streams adopting studies carried out by Rôdas, (2013) and Garcia (2013)  We used 13 years of water quality data to estimate the monthly geometric mean of concentrations based on a recent approach suggested by Isles (2020).
Daily flows measured in the streams were made available between December 2013 and December 2017 by two fluviometric monitoring stations located in the Sorocabuçu and Sorocamirim streams (62472800, 62473200). We calculated a mass balance to determine the reservoir inflow concentrations ( ) using Eq 1:  Rôdas (2013) and Garcia (2013), and monitoring stations operated by CETESB and ANA.
Particulate organic carbon was calculated as the difference between the total organic carbon and the dissolved organic carbon. Organic nitrogen was calculated as the difference between the Total Kjeldahl nitrogen (TKN) and ammonium (NH4) concentrations. We considered that the organic nitrogen concentrations would be higher in the particulate organic nitrogen (70%) than in the dissolved organic nitrogen (30%), based on the stream's features.
We applied specific ratios to individual phosphorus forms based on Garcia The input chlorophyll-a (Chla) concentrations followed the same mass balance equation using a conversion factor of 50mgC.mgChla-1.

Biogeochemical modeling
The General Lake Model (GLM) is a one-dimensional hydrodynamic model that calculates vertical profiles of water temperature, salinity and density by representing upstream and downstream water flows, mixing, heating and surface cooling (Hipsey et al., 2017). The model can be coupled with the Aquatic Ecodynamics Modeling Library (AED2).
The AED2 comprises modules and algorithms that simulate water quality, aquatic biogeochemistry, and phytoplankton and zooplankton dynamics in lakes, reservoirs and estuaries (Hipsey et al., 2013). The modules allow the simulation of carbon, nitrogen, phosphorus, dissolved oxygen and silica cycles, as well as organic matter, phytoplankton and zooplankton through mass balance equations.
In the scope of the present chapter, the mass balance equations of dissolved oxygen (DO) and nitrate (NO3) concentration calculations are detailed below: Where: 2 =atmospheric O2 exchange, 2 =sediment O2 demand, consumption by zooplankton respiration.
Where: ̂ = and is the thickness of the z th layer/cell Where: ]=uptake from the phytoplankton community

Calibration and validation
The model parameters were calibrated following the bottom-up principle: firstly, the water level and water temperature were previously calibrated as described in Barbosa et al. (2021), and then the parameters regarding the concentrations of DO and NO3. Due to the low water quality data availability compared to the water level and water temperature historical time-series, the chosen calibration and validation period were shorter for the present study. We A global sensitivity analysis based on the Morris Method (Morris, 1991) was performed to identify the most sensitive parameters for the predictions of DO and NO3. After that, an automatic calibration was performed using the derivativefree, optimization algorithm (CMA-ES; Hansen, 2016) with 100 iterations aiming to reduce the root mean square error (RMSE) followed by a manual calibration to ensure that the model was not reproducing unreal biogeochemical parameter combinations. A similar sensitivity analysis and calibration approach was conducted by Ladwig et al. (2020).
The model parameters were manually changed aiming to sequentially optimize goodness-of-fit (GOF) metrics focusing on reproducing DO and NO3 concentrations of the dry and wet periods. As the observed NO3 concentrations in the reservoir did not show a significant temporal fluctuation pattern during the calibration and validation periods, we focused on representing mainly the longterm median concentration. As the purpose of this study was to simulate future scenarios and represent the likely alterations and quantify them based on a baseline scenario, we did not focus on capturing NO3 peaks during the calibration process.
We calculated five GOF metrics (RMSE, the Pearson correlation coefficient (r), mean absolute error (MAE), the Nash-Sutcliffe model efficiency coefficient (NSE) and the Kling-Gupta efficiency (KGE)) to compare model outputs regarding the profile, surface and epilimnion concentrations and measured data using the hydroGOF package for R. (Zambrano-Bigiarini, 2017).
The target of calibration was to maximize the r and reduce the RMSE values for each variable. A similar approach was performed by Fenocchi et al. (2019).

Basin distributed load modeling
We adopted the methodology developed by Anjinho et al. (2021) based on the export coefficient modelling approach (Johnes, 1996)  pixel. We adopted 13.5 m³.s as the long-term mean daily inflow in the reservoir  to simulate the accumulated long-term mean annual flow per pixel for each upstream. We used the regionalization method based on a basin yield that assumes the existence of a proportional linear relationship between drainage area and streamflow to simulate distributed inflow in the basin.
Thus, we divided the long-term mean daily inflow by the total number of pixels in the basin (m3 s−1 pixel−1) and the flow accumulation algorithm was used to determine the model of mean annual accumulated inflow. kg.km -2 .y -1 to kg.pixel -1 .y -1 and then the flow accumulation algorithm was used to generate the accumulated nutrient load model.

Land use scenarios for the 2050s
Three land use scenarios were considered to evaluate likely future impacts on the DO and NO3 concentrations in the Itupararanga reservoir.
The potential transition scenario (MOLUSCE) was used as a first scenario and also a baseline for future modification actions in the basin. The other two scenarios were formulated based on the simulated LULC for 2059 considering the increase in preservation areas, focusing on the restoration of permanent preservation areas and reducing agricultural uses (Green Scenario) and increased economic development of the basin focusing on agriculture, pasture, soy and urbanization (ED scenario) ( Table 2). The inflow mean TP and TN concentrations were accounted for in each scenario and those concentrations were compared to the baseline period. Thus, we altered the NO3 inflow and FRP concentrations by the increased or decreased proportions based on the baseline values in the biogeochemical model. This modelling approach aimed to take into account the influence of the LULC changes in the water quality of the Itupararanga reservoir.

Biogeochemical simulations
The input time-series concentrations show that overall, the pH values and DO concentrations represented peaks at the beginning of the dry season (June).
On the other hand, peaks of TN concentrations were at the end of the dry season (September). The TOC and TP concentrations followed the expectations with higher values in the wet season due to the watershed runoff increase ( Figure S1, Supplementary material).
The biogeochemical sensitive parameters and their respective values calibrated in the present study are compared with the reference values and shown in Table S1. Overall, the biogeochemical modeling showed good fit criteria compared with the observations for both simulations of the DO and NO3 concentrations ( Table 3). The DO simulation was able to catch the small changes between wet and dry season, especially in the calibration which has a longer period of simulation compared to validation (Figure 3).   The validation results of the transition potential for the Alto Sorocaba basin in 2019 using 1999 and 2009 LULC data according to each land use/land cover are presented in Table 5 based on percentage areas in relation to the total area of the Alto Sorocaba basin.  (baseline) is shown in Figure 4 and the built scenarios related to the baseline are shown in Figure 5.

Changes in DO and NO3 concentrations under three land use scenarios
The DO concentrations from the three scenarios slightly changed compared to the baseline (Table 6, Figure 6a). As expected, the mean DO concentration decreased in the MOLUSCE and ED scenarios due to an increase in O2 consumption by nitrification. On the other hand, changes in inflow NO3 and PO4 concentrations led to a greater influence on the NO3 simulations in the reservoir, which may be associated with the greater impact of allochthonous loads in the reservoir due to changes in land use (Figure 6b). The simulated Chla concentrations for each scenario were shown to be influenced by such changes of the inflow dissolved nitrogen and phosphorus concentrations (Figure 6c).

Performance of models
The biogeochemical modeling was able to represent the average DO and c. The distributed modeling validations of TN and TP concentrations along the Alto Sorocaba Basin based on the nutrient loads exported from the catchment (Johnes, 1996) can highlight the potentialities of this GIS approach to assess the LULC impacts on the streams and the Itupararanga reservoir. The simulated TP and TN concentrations in the watercourses showed good fit criteria compared with the reference values from Moriasi et al. (2007). Another recent study applied a catchment-scale nutrient model with a similar modeling fit compared to the results shown in Table 4 (Messina et al., 2020).
Despite the good results generated by this modeling approach in the present and previous studies (Anjinho et al., 2021;Lima et al., 2016), it represents the dynamics of nutrients in rural basins more effectively than in urban ones due to a poor performance in the simulation of nutrient concentrations from point source pollution. Since the Alto Sorocaba Basin has less than 10% of the urban area in the total catchment area, we did not consider point sources as sewage treatment plants in the basin for our simulations. Another limitation of this approach is considering only the conservative nutrient transport that does not take into account the temporal and spatial transformations of TN and TP concentrations in watercourses (Lima et al., 2016).

LULC changes in the last two decades and projections for 2050s
We have analyzed LULC changes in the Alto Sorocaba basin comparing the last two decades (1999 and 2019) based on data from the MapBiomas project (2020). It can be observed that the pasture areas reduced 50% of their area in that period, but the agricultural areas increased significantly focusing on sugarcane, which grew three times, and soybeans, which grew five times its area, in addition to forest plantation areas (which grew four times their previous percentage of area). The observed LULC changes were similar to the spatial pattern already highlighted for Brazil in the past decades (Miccolis et al., 2014).
The APA Itupararanga has favored environmental conservation in the Alto has been widely used to assess and predict LULC changes from small (Satya et al., 2020) to large areas (Fernandes et al., 2020). Taniwaki et al. (2013) suggested that exposed bare soil, agriculture without the protection of riparian zones and urbanization were the most significant drivers of water quality degradation and reflected major damage to the Itupararanga reservoir in 2010. The low availability of NO3 concentrations and the phosphorus limitation in the Itupararanga reservoir has proven to have significant effects on cyanobacteria biomass, especially R. raciborskii, and toxin levels measured in the reservoir surface waters (Machado et al., 2021, under review (Curtarelli et al., 2015;Lins et al., 2018;Ma et al., 2016).
We assess water quality changes focusing on reservoir DO and NO3   (BRASIL, 2005), based on the kernel density estimation (KDE) (Figure 6c), the average chlorophyll-a concentrations would be above the limit (<30ug.l) in the ED scenario which may indicate alterations in the trophic dynamics of the reservoir.
Higher Chla concentrations were expected in such a scenario compared to the others, mainly because phosphorus is the final limiting nutrient for primary production in (sub)tropical freshwaters (Quinlan et al., 2020). However, the high potential for NH4 and NO3 uptake by cyanobacteria throughout different periods of the year has been reported in the Itupararanga reservoir (Cunha et al., 2017), mainly driven by the N availability and high water temperatures. It is also important to consider warming air temperatures and other climate change projections to assess the likely consequences for reservoir water quality, as suggested by Barbosa et al., (2021).

CONCLUSION
The present study aimed to simulate future LULC scenarios in the Itupararanga reservoir based on likely conditions of LULC in 2059 and two hypothetical conditions considering an increase in preserved areas (green) and another accounting for the increase in economic development activities focusing on agricultural and urbanized areas.
Thus, we calibrated and validated a biogeochemical model focusing on dissolved oxygen and nitrate concentrations to simulate their responses from the LULC scenarios. Even if DO and NO3 concentrations did not increase above the limits allowed by Brazilian legislation for all scenarios, Chla concentrations would increase significantly, which could pose a serious threat to the use of the Ituparanga reservoir as a source of water supply. The coupling modeling was also able to represent the importance of forest formation lands in the reservoir water quality improvement (green scenario). We did not assess the management characteristics in relation to exposed soil and agriculture without riparian protection in detail in the present study, however, we suggest that these characteristics should be considered in future studies on LULC scenarios in the Itupararanga reservoir.
The modeling framework proposed based on the coupling watershed load and biogeochemical lake model was able to represent the responses of water quality changes in the reservoir due to the LULC changes. This novel modeling framework can be a valuable tool to guide water resources management considering future pressure on freshwater due to population growth and intensive agricultural practices for human consumption.
The basin distributed load approach is a GIS implementation of commonly used simple nutrient loading equations and can be run in any GIS software.
Detailed instructions are provided in Anjinho et al., 2021. Zambrano-Bigiarini, M., 2017. Package 'hydroGOF'. In: Goodness-of-fit Functions for Comparison of Simulated and Observed. Figure S1: Input time-series of pH and DO, TOC, TP and TN concentrations used in the biogeochemical simulations.