Impact of large reservoirs on simulated discharges of Brazilian rivers

Tens of thousands of dams were built around the world to reduce flood risks, produce energy, and maximize benefits of limited freshwater resources. In Brazil, the main and largest reservoirs are related to hydropower plants. Improving the understanding of reservoir dynamics is important not only to evaluate their impact in the flow regime of Brazilian rivers, but also to simulate the combined effect of constructing new dams and potential alterations under future climatic conditions. Here, we analyze how an ideal representation of reservoirs in terms of forced discharge would improve a previously calibrated hydrological model under the Brazilian domain. We forced the continental-scale version of the MGB model on observed reservoir outflows from 109 hydropower dams, which are part of the Brazilian National Interconnected System controlled by the National Electrical System Operator. Model simulated flows were replaced by the reservoir outflows in all dam locations and were compared to the original discharge in downstream gauges. The forced discharge simulation presented a mean improvement for Kling-Gupta Efficiency of 21%, when compared to the original model (naturalized flow). This analysis is a preliminary step towards an explicit representation of the reservoirs in the model, what will be conducted in a future study.


INTRODUCTION
The hydrological cycle of the Earth is changing due to anthropogenic alterations in water resources management and land uses and cover. In one hundred years, the population has quadrupled to 7 billion people, and while half of them live in urban areas, rapid urbanization is occurring in many developed and developing regions of the world (Klein Goldewijk et al., 2011). For food, the irrigated agricultural area increased six times in the last century (Siebert et al., 2015) and the water demand continues to rise, driving an ever-increasing agricultural demand (>70% of all fresh water). In Brazil, irrigation represents about 50% of the total estimated abstractions and 70% of the total estimated water consumption (Agência Nacional de Águas, 2018).
In order to maximize the benefits of limited freshwater resources and mitigate flood risks, tens of thousands of artificial dams have been built in the world's major river systems (Nilsson et al., 2005;Lehner et al., 2011). Hydropower generation is the main use for the largest reservoirs in Brazil, corresponding to 65% of the country's installed capacity (Empresa de Pesquisa Energética, 2018). Due to the severe droughts that hit Brazil's main basins between 2012 and 2015, the average share of hydroelectric power generation declined to 76% due to nation scale depletion of reservoirs, with great use of thermal plants and large reservoir depletion. Just before, between 2000 and 2012, artificial dams yielded an average of 91% of the power generation (Zambon, 2015).
Despite of the national relevance of large reservoirs for energy production and reduction of flood risk, it is fundamental to understand how these major structural interventions change natural hydrological regimes. In general, reservoirs store water and release it downstream according to the current demand. Regulating water flow results in attenuated hydrographs, higher evaporation due to the increase in water surface area, and also major changes in sediment and nutrient transport (Ligon et al., 1995). The existence of dams also leads to fragmentation of rivers, reducing the natural connectivity within and between fluvial systems (Tischendorf & Fahrig, 2000;Moilanen & Hanski, 2001), among other impacts.
An interesting way to assess the impact of dams on river dynamics is by making use of conceptual and physically-based hydrological models. In this context, reservoirs can be simulated explicitly by incorporating them in the model structure (e.g., solving the water balance for storage with stage-area-volume curves and operation rules) or not explicitly by using a simple substitution form, in which the simulated flows are replaced by reservoir operation data at the corresponding dam locations (i.e., the hydrological model is forced with observed reservoir outflows). In Brazil, reservoirs have been simulated in hydrological models either by explicit or substitution approaches (e.g. Fleischmann et al., 2019bFleischmann et al., , 2019c, but the existing studies usually have a focus in specific river basins or particular dams. A comprehensive analysis regarding the impact of reservoirs at the national scale, on the other hand, has not been carried out so far. Advances in modeling techniques, computer capabilities and big data has taken place in the last years (Siqueira et al., 2018), and assessing reservoir impacts over large spatial domains has become even more a reality. For instance, in Global Hydrology and Water Resources models (GHWRM), which aim to assess the anthropogenic impacts over the entire globe (see Bierkens, 2015), numerous generic algorithms have been developed to explicitly represent the influence of lakes and reservoirs on river flows (Döll et al., 2003;Haddeland et al., 2006;Hanasaki et al., 2006;Wada et al., 2014;Zajac et al., 2017;Shin et al. 2018). However, generic operation rules used by GHWRMs are rather simplistic and often rely on global datasets (e.g., Lehner et al., 2011), which may hinder the representation of complex reservoir systems such as the Brazilian National Interconnected System (SIN). Exploring regionally available dam outflow data using a hydrological model with national coverage would be an important step to overcome such limitations, helping further development for robust operation rules and giving insights on the potential benefits of including reservoirs in the model structure. The main Brazilian dataset is the Reservoir Monitoring System (Sistema de Acompanhamento de Reservatórios) from the National Water Agency (ANA, in the Portuguese acronym), which provides daily time series of information as water level, percentage of active storage, inflows and outflows.
MGB (Collischonn, 2001;Collischonn et al., 2007;Paiva et al., 2013) is a large scale hydrological model which was already applied to the whole South America in a multi-basin manner with satisfactory performance (Siqueira et al., 2018). Although some studies have explored the explicit simulation of reservoirs in MGB (Larentis et al., 2008;Collischonn et al., 2011;Fleischmann et al., 2015Fleischmann et al., , 2019bFleischmann et al., , 2019c, reservoir modeling is not yet consolidated in MGB. Then, the aim of this work is to evaluate how an ideal representation of SIN reservoirs in terms of forced discharge would improve a continental-scale, previously calibrated hydrological model.
The paper is divided as follows: firstly in the Material and Methods section, the model is described, together with reservoir data and performance evaluation criteria. Then results are shown for an exploratory analysis of the reservoirs behavior, followed by the model with forced discharges, and finally for the differences in discharges between natural and forced simulation scenarios.

Model description
MGB is a conceptual, semi-distributed hydrological model designed to simulate land surface hydrological processes and flow propagation in tropical basins (Collischonn et al., 2007). It divides the basin into unit-catchments to capture the spatial variability of both meteorological and physical characteristics and further into hydrological response units (HRUs) to represent land use and soil heterogeneities. Water and energy budget are computed at the HRU level considering the soil as a bucket model. Evapotranspiration (from soil, plant and canopy) is calculated using the Penman-Monteith equation. Surface, subsurface and groundwater runoff generated within the unit-catchment are propagated to the stream network using linear reservoirs to represent the catchment damping effect. Flow routing in river channels can be computed using different approaches, ranging from simple methods such as the Muskingum-Cunge (Collischonn et al., 2007) to more complex ones like the full 1D St. Venant equations (Paiva et al., 2013). Passaia et al.

3/9
The continental MGB model version for South America (MGB-SA) (Siqueira et al., 2018) is used in this study. It uses the inertial equation of propagation in rivers described in Pontes et al., (2017), with sufficient physical basis to represent river hydrodynamics and floodplains effects. This version has 33,749 unit-catchment with river reaches at every 15 km and adopts the Multi-Source Weighted Ensemble Precipitation (MSWEP) as the forcing data, a product that combines estimates of rainfall by satellite, reanalysis and observed precipitation . The model was manually calibrated for the period between 1990 and 2010 with >600 discharge gauge stations (including reservoir naturalized flows from ONS, i.e., without the effect of regulation), and was validated using multiple remote sensing databases. Performances for daily flow regarding the Kling-Gupta and Nash-Sutcliffe efficiencies were > 0.6, respectively, for 70% and 55% of the gauges (Siqueira et al., 2018). Here we used the same simulation period as the original model version, from 01/01/1990 to 31/12/2009.

Reservoir data and substitution of flows in the model
To evaluate the impact of reservoir regulation on water discharges, the model was forced with observed outflows from reservoirs controlled by the ONS. For that, simulated discharges at unit-catchments corresponding to reservoir locations were replaced by the observed outflows at each model time step (daily), and this simulation was called the F simulation. Reservoirs that meet one or more of the following criteria were not considered in this study: a) built near the end simulation date of the model; b) located in the same unit-catchment and upstream of another reservoir, since in the approach used here only one reservoir per unit-catchment is possible (the most downstream reservoir was simulated); c) outside the simulated spatial domain; and d) with insufficient flow data. Thus, it was included discharge data from 109 reservoirs (out of an original dataset of 150 reservoirs), obtained from the Reservoir Monitoring System (SAR, in the Portuguese acronym). Figure 1 shows the name and location of all reservoirs used in this study. Only the impact on the flows was analyzed.
We carried out an exploratory analysis of the reservoir SAR data by a visual inspection, and also computed the correlation between the outflows and other variables: inflows, storage (water level and volume) and the day of the year (Julian day). For the latter case, the correlation was computed in terms of circular statistics to avoid the discontinuity of days 365 and 1.

Evaluation of model performance
The F simulation and naturalized scenarios (i.e., the naturalized flow simulation without reservoirs, called here N simulation) were compared with in situ observed flow data obtained from the National Water Agency (ANA) monitoring network.
The Kling-Gupta efficiency indexes (KGE, Gupta et al., 2009) and Nash-Sutcliffe (NSE, Nash & Sutcliffe, 1970), correlation (r) and volume error (bias) were also used. Another metric used by Zajac et al. (2017) is the skill score (IMD, in the Portuguese acronym). A positive IMD means that the performance of the model with reservoirs has improved, while negative values indicate the opposite trend. IMD can be applied to any performance index; an example is given in Equation 1 by computing IMD for the KGE.
The skill score for the KGE (IMD KGE ) is a function of the KGE for the reservoirs simulation (KGE F ), the KGE for the naturalized flows simulation (KGE N ), and the optimum KGE (KGE opt = 1).
Changes in flows were analyzed with the following metrics: normalized root mean square error (nRMSE -i.e., RMSE divided by average discharge), correlation (r), standard deviation, and absolute differences in Q 10 and Q 90 , between the N and F simulation. The average absolute percentage difference between N and F simulations is computed annually for the main basins of South America (Amazon, Tocantins, São Francisco and Paraná).

Exploratory analysis of the reservoirs' behavior
The average behavior of some variables throughout the year was investigated. Three dams were chosen to cover different regions and volume capacity (Figure 1). Based on Figure 2b it is possible to infer that, for the Porto Estrela reservoir, there is a period with lower flow release between May and mid-October. For all analyzed reservoirs, the outflow can vary expressively for The labeled reservoirs are those whose results will be discussed in the text. the same day, in different years. The dispersion is visually large, for both outflow and inflow (Figure 2a), due to the climate variability. The correlations were in average 0.2, indicating a certain relation between the day of the year and the outflow.
One can conclude the same from the outflow to the inflow. The largest difference is in the Sobradinho reservoir correlation. Its correlation for inflow is 0.79; and for outflow, 0.16, what shows its great storage capacity and the impact on downstream flow regime. In seasonal basins, a greater correlation between inflow and season is expected.
Regarding the relationship between daily water elevation and outflow (Figure 2c), it is hard to deduce any operating rule based on its behavior. That was expected because ONS uses a complex operation rule, not based only in the reservoir's water elevation. Moreover, the analysis performed here considered only daily data; a future assessment based on monthly data could lead to different conclusions. In the figure, it is observed that the reservoirs of Porto Estrela and Sobradinho have their highest discharge values in summer and fall, while 14 de Julho, a run-of-the-river plant, in winter and spring.
The analysis of the relation between outflow and possible driving factors (e.g. storage, inflow, time of year) can be studied through the correlation coefficient (Figure 3). It was done for 109 reservoirs. The variable that most explains the outflow is the inflow. The mean correlation is 0.79. The second most relevant factor is storage, represented here by volume and water elevation, with mean correlations of 0.31 and 0.22, respectively. The circular correlation between outflow and day of the year presented a mean correlation of 0.19, being the least important factor. Figure 4 presents a spatial assessment of the IMD metric for KGE index. In most gauge stations (≅ 74%) the performance improved. A poorer performance was observed along the southeast Figure 2. (a, b) daily inflow/outflow during the year, for three Brazilian hydroelectric reservoirs after their constructions. The average for all the SAR reservoir's data is shown in red, all observations in blue; in the title, the reservoir's name and the correlation between the inflow/outflow and the day of the year. Some points are not shown to better visualization of the pattern; (c) outflow versus storage (water elevation) of three Brazilian hydroelectric reservoirs, comparing the years seasons.

Forced Discharge Simulation
Passaia et al.

5/9
and extreme west coast of Brazil (Paraíba do Sul and Alto Paraguai basins). This may be caused by uncertainties related to input data in all the domain. Furthermore, since the MGB-SA was calibrated for the situation without reservoirs and, although care has been taken to use natural observed discharge and gauge stations without reservoirs influence, the calibrated parameters (such as the soil storage capacity parameter or the relationship between storage and saturation in the soil parameter) may have been adjusted in the natural scenario with the effect of the reservoirs on the outflow. That is, the effect of the reservoirs in the simulation was possibly compensated during model parameterization. After forcing the model with a reservoir outflow, it would be necessary to recalibrate the model to remove this compensation. This is the main reason we may have lower results in the forced discharge scenario at some gauge stations.
Due to MGB-SA inability to perform well in the Brazilian Northeast (semi-arid climate, Kottek et al., 2006), there was a considerable improvement in the results downstream of São Francisco river, mainly due to the consideration of Sobradinho and Luiz Gonzaga reservoirs, close to the river mouth. The Três Marias reservoir discharge also has an important role there. Most of the discharge in São Francisco River is generated in these south headwaters, since the downstream reaches are located in a more arid region. The impact of Brazilian reservoirs on discharge can be seen even outside Brazil, affecting other countries such as Paraguay and Argentina, located in downstream parts of the Paraguay and Paraná rivers (Bonetto et al., 1987).
There was no clear relationship between IMD and reservoirs volumes ( Figure 5). For volumes larger than 1000 hm 3 there seems to be a trend towards IMD improvement with increasing  volume. The F simulation showed improvement in 80% of the stations and, of these, 40% improved more than 20% (IMD> 0.2). Hydrographs downstream of three selected reservoirs are presented in Figure 6 to compare the F and N simulations (see Figure 1 for dams' location). The reservoir of Irapé (Jequitinhonha River) presented significant improvement in the hydrograph, whereas, for the others, the model with naturalized flows was satisfactory. The 14 de Julho dam in Taquari-Antas River is a run-of-the-river plant, which almost does not alter the downstream flow regime. Great improvement is observed in the base flow downstream of Manso reservoir in the Cuiabá River, which was underestimated in the naturalized flows model.
The calculated metrics for discharge gauge stations downstream of the 109 dams were summarized by the 10th percentile, median and 90th percentile (Table 1). F scenario had a performance superior than scenario N. The metric with the greatest improvement was the logarithm of the NSE (due to more uncertainties in the model base flow), followed by the NSE, KGE and r. The KGE showed an improvement of 21% for the F simulation.
Each metric depicts different flow signatures. The correlation reflects the seasonality, but is unable to evaluate the average flow.
A NSE equal to zero shows that the simulations are as accurate as the mean of the observed data, and negative NSE means that the mean is a better estimate than the model. To further improve the assessment of hydrological models, the KGE was created (Gupta et al., 2009), which takes into account the correlation, bias and variability. This way, it is ensured that bias and variability are not correlated.

Differences in discharges between natural and forced simulation scenarios
Changes in daily discharges between simulations F and N were evaluated through the following metrics: nRMSE, and relative mean absolute differences in Q10 and Q90 (Figure 7a, b and c,   respectively). The difference between both scenarios must be analyzed considering that MGB-SA estimates present relevant biases in certain basins, which is registered in the supplementary material of Siqueira et al. (2018) and is summarized in Figure 7d for the 109 dams. In this case, we computed the bias between scenario N and the observed outflow for each dam. For a future study, the impact of reservoirs on river discharges may be analyzed by firstly performing a bias correction on the model naturalized flows. The nRMSE calculated exceeded 80% in the lower part of the São Francisco River, while MGB-SA bias is lower than 50% for most gauges there (Figure 7d), indicating that although the model has a large bias in such semi-arid regions, important flow regulation is performed by the reservoirs. The Paraná River presents nRMSE higher than 60% in its upstream portion, nRMSE between 20 and 40% until the Brazilian border, and less than 20% in its downstream reaches. Compared to the Paraná River MGB-SA bias (smaller than 10% for most dams, Figure 7d), nRMSE indicates that considerable portions of the basin upstream from the Brazilian border are under flow regulation. The Amazon basin is little affected by reservoirs in terms of flow.
In general, the difference in Q10 was less than 20%, while the difference in Q90 was greater for the Paraná River and its tributaries (20-40%). However, we cannot say that this difference is due to the reservoirs' presence, because the MGB-SA flow bias is, on average, greater than 20%.

CONCLUSIONS
The aim of this study was to evaluate how an ideal representation of reservoirs in terms of forced discharge would improve a continental-scale and previously calibrated hydrological model. For this purpose, we considered outflows of 109 SIN reservoirs in the MGB-SA model (Siqueira et al., 2018).
Through an exploratory analysis, it is identified how observed reservoirs outflows are related to inflows and water elevation. As most reservoirs are run-of-the-river (64 out of 109), there is a strong relation between inflow and outflow, however, no expressive correlation is observed between outflow and water storage. Thus, it is not trivial to build a general operation rule for all 109 SIN reservoirs.
The simulation with forced discharge (F simulation) presented better performance than the MGB-SA with naturalized flows (N simulation), showing a 21% improvement in KGE. Large rivers as São Francisco and Paraná were shown to be largely regulated, even considering the MGB-SA bias. In a next step, we intend to explicitly represent the reservoirs in the model and test operating rules. It is possible to develop a model of water uses, and then use the output schemes developed to meet the demand of water downstream. In addition, level-volume relationships provided by the ONS and the new ANA water use products can be used (Agência Nacional de Águas, 2018).

ACKNOWLEDGMENTS
The present work was carried out with the support of CNPq (National Council of Scientific and Technological Development /Brazil), under the grant numbers 132663/2017-1 and 141161/2017-5. The authors also acknowledges IPH and UFRGS for their support, and the anonymous reviewer's suggestions that significantly improve the text.