Large-Scale Hydrological Modelling of the Upper Paran á River Basin

: The Upper Paran á River Basin (UPRB) has undergone many rapid land use changes in recent decades, due to accelerating population growth. Thus, the prediction of water resources has crucial importance in improving planning and sustainable management. This paper presents a large-scale hydrological modelling of the UPRB, using the Soil and Water Assessment Tool (SWAT) model. The model was calibrated and validated for 78 outlets, over a 32-year simulation period between 1984 and 2015. The results and the comparison between observed and simulated values showed that after the calibration process, most of the outlets performed to a satisfactory level or better in all objective functions analyzed with 86%, 92%, 76%, 88%, and 74% for Percent bias, Coe ﬃ cient of determination, Nash-Sutcli ﬀ e e ﬃ ciency, Kling-Gupta e ﬃ ciency, and the Ratio of Standard deviation of observations to root mean square error, respectively. The model output provided in this work could be used in further simulations, such as the evaluation of the impacts of land use change or climate change on river ﬂows of the Upper Paran á Basin. ﬂow e ﬀ ective chosen.


Introduction
Hydrological models have been used worldwide as a powerful tool for water resources research and management. Many studies have focused on modelling the hydrology of areas that have experienced an increase in the frequency of drought and flood events. Over recent decades, hydrological modelling has contributed to improving the conservation and sustainable use of natural resources, especially through research activities dedicated to mitigating climate change [1,2], land use changes [3,4], and sources of water pollution [5]. However, most such studies are dedicated to small-to medium-sized basins, which produce some difficulty in generalizing the conclusions to large-scale basins. Notably, collecting and organizing a good set of data that describes the physical properties of a small river basins well, is considerably easier than doing the same for a large river basins. The challenge in preparing input data, with high spatial and temporal resolutions, is another factor that prevents hydrological modelling studies from focusing on large-scale river basins. Therefore, only a few studies have been performed for large-scale basins [6].
The Upper Paraná River Basin (UPRB), located in central-southern Brazil, is one of the largest and most socio-economically important basins in South America. It plays a significant role in the Brazilian economy and development, greatly contributing to economic sectors, such as agriculture, livestock, energy, and urban and industrial water supply. In particular, this watershed houses 87 hydropower Before reaching the border between Brazil and Paraguay, the Paraná river receives large and socio-economically important tributaries on the east side of the basin, such as the Tietê and the Paranapanema rivers, in São Paulo state. In addition, the west side of the Paraná river crosses the Maracaju mountain range, in Mato Grosso do Sul, which acts as a natural barrier separating the Pantanal wetlands and leads to the formation of many rivers, shorter than those in the east side of the basin.
The UPRB has a unique geographical profile, with a number of hydropower plants close to the largest urban and industrial areas, which are large consumers of electricity. The basin has an estimated population of over 65 million inhabitants, of whom more than 93% live in urban areas [14]. According to the Brazilian National Water Agency (ANA) [15], this region has the highest demand for water resources in Brazil, equivalent to 736 m 3 s −1 , mostly used for agricultural (42%) and industrial (27%) activities.
The UPRB is embedded within the center-east portion of South America, with an approximately oval shape, and with the major axis in the north-south direction. The basin is characterized by different morphologies that range from Atlantic Plateau (elevation higher than 2000 meters) to the Paraná River Valley (between 350 and 100 meters). It is a sedimentary and igneous basin, with the volcanic rocks of the Serra Geral formation overlaid by sedimentary rocks, mostly located in the central and western region of the basin [16][17][18]. Sedimentary areas are also found in the contours of the basin, in their higher hills. This type of formation, combined with volcanic rocks, predominates in most of the tributaries and progresses up to near the main course of the UPRB. Most areas of basaltic rocks are formed by high fertility soils. Until a century ago, such areas were covered by a dense forest, predominantly with medium-to-large tree. This forest cover was almost completely removed within the basin and the exposed land was replaced by intensive agricultural exploitation. During the 20th century, rapid population growth and pronounced urbanization have led to a significant land use change in the UPRB. For instance, Paraná and São Paulo states, located in the east of the basin, have lost more than 70% of their primitive forest, while the original vegetation in the western part of the basin, was maintained until the 1970s, when the development of agro-business increased. Deforestation occurred for different objectives, but in most cases, forests were replaced by agriculture and pasture [11].
The different land covers and the intense internal dynamics of the land uses may have affected the regional hydrology in different ways, since some areas of the basin have increased, while other areas have decreased their stream flows [12,13]. Therefore, studies on the simulation of water resources in the UPRB have great importance in offering subsidies for managers and policy makers. A better understanding of the UPRB hydrology could improve the planning and sustainable management of the wide range of water uses in the basin.
Considering the importance of the UPRB and its significant changes in stream flow, since the mid-1900s, the primary goal of this work was to use the Soil and Water Assessment Tool (SWAT) model to estimate the discharge in monthly time steps at a highest spatial resolution allowed by the simulation system. To achieve this goal, this work pursued the following main objectives: (a) set up the SWAT model with the most appropriate dataset available; (b) calibrate and validate the main outlets of UPRB, including uncertainty assessment; (c) evaluate the performance of the model using several objective functions; and (d) address the spatial and temporal analysis of discharge over the UPRB.
The results of this work address the variability of the discharge at the basin as a whole in the spatial and temporal dimensions. Additionally, the approaches and strategies used for calibration might serve as standards for future simulations of large river basins. Furthermore, this work creates a basis for future studies on the UPRB to assess the potential impacts on hydrological processes and water quality, allowing the simulation of diverse scenarios such as climate and land use changes.
The remainder of this paper is structured as follows. First, it gives a brief description of the study area, the hydrological model used for simulation, the input data, and the setup performed to build a project for UPRB. Then, the strategies for calibration, sensitivity and uncertainty assessment are described in the following three sub-sections: SUFI-2 and parameter calibration, objective function, and modelling protocol. Finally, the results are presented and discussed and are followed by conclusions.

Study Area
The UPRB is located in the central-southern region of Brazil ( Figure 1), with an area of 900,480 km 2 , and drains rivers in six Brazilian states: São Paulo (23.5%), Paraná (20.4%), Mato Grosso do Sul (18.9%), Minas Gerais (17.6%), Goiás (15.7%), Santa Catarina (1.2%), and the Federal District (0.4%), as well as a small portion of Paraguay (2.3%). The Paraná river is the second largest river in South America. From the confluence of the Paranaiba and Grande rivers, the Paraná river flows southward for 738 km until it reaches the border between Brazil, Argentina, and Paraguay.
Before reaching the border between Brazil and Paraguay, the Paraná river receives large and socio-economically important tributaries on the east side of the basin, such as the Tietê and the Paranapanema rivers, in São Paulo state. In addition, the west side of the Paraná river crosses the Maracaju mountain range, in Mato Grosso do Sul, which acts as a natural barrier separating the Pantanal wetlands and leads to the formation of many rivers, shorter than those in the east side of the basin.
The UPRB has a unique geographical profile, with a number of hydropower plants close to the largest urban and industrial areas, which are large consumers of electricity. The basin has an estimated population of over 65 million inhabitants, of whom more than 93% live in urban areas [14]. According to the Brazilian National Water Agency (ANA) [15], this region has the highest demand for water resources in Brazil, equivalent to 736 m 3 s −1 , mostly used for agricultural (42%) and industrial (27%) activities.
The UPRB is embedded within the center-east portion of South America, with an approximately oval shape, and with the major axis in the north-south direction. The basin is characterized by different morphologies that range from Atlantic Plateau (elevation higher than 2000 m) to the Paraná River Valley (between 350 and 100 m). It is a sedimentary and igneous basin, with the volcanic rocks of the Serra Geral formation overlaid by sedimentary rocks, mostly located in the central and western region of the basin [16][17][18]. Sedimentary areas are also found in the contours of the basin, in their higher hills. This type of formation, combined with volcanic rocks, predominates in most of the tributaries and progresses up to near the main course of the UPRB. Most areas of basaltic rocks are formed by high fertility soils. Until a century ago, such areas were covered by a dense forest, predominantly with medium-to-large tree. This forest cover was almost completely removed within the basin and the exposed land was replaced by intensive agricultural exploitation.
The basin has great spatial variability in its precipitation pattern as it covers different climatic areas. Climatologically, the northern part of the UPRB is influenced by the South American Monsoon System [19,20], characterized by wet summers and dry winters. Precipitation during the summer may exceed 800 mm, while during winter it can be as low as 30 mm. In the southern parts of the basin, precipitation is spread out throughout the year and is associated with the Intertropical Convergence Zone (ITCZ), cold fronts, and the Mesoscale Convective Complex (MCC), mainly during the spring and summer [21,22].

SWAT Model
The hydrological simulations of UPRB were performed using the 2012 version of the Soil and Water Assessment Tool (SWAT) model with an ArcGIS interface. SWAT is an open source, semi-distributed, and physical model developed by the Agricultural Research Service of the United States Department of Agriculture (ARS-USDA). This model can be used to design analyses related to physical processes, both in small and large watersheds, and can be executed in a continuous simulation in monthly or daily time steps. It is widely used to assess impacts on hydrological processes, water quality (e.g., transport of nutrients and pesticides), as well as climate and land use change scenarios [23][24][25]. Based on the topography, a basin is discretized into sub-basins, which are connected by a stream network. To assess the differences in land cover and the heterogeneous soil in a watershed, each sub-basin is further discretized into hydrologic response units (HRUs), according to unique combinations of land use, soil type, and slopes. For each HRU, simulated hydrological processes, such as surface runoff and evapotranspiration, are generated separately, and then routed through the river network to the outlet of the basin. For further details on the SWAT model, the reader is referred to Neitsch et al. [26].

Data Description and Model Set Up
Different input data are required to build a hydrological project with SWAT, including meteorological, hydrologic, and physical variables. The data used in this work was prepared for the whole simulation period (which includes the warming up, calibration, and validation periods), from January 1979 to December 2015. The first five years (1979)(1980)(1981)(1982)(1983) were used for warming up, the following 21 years  constituted the calibration period, and the last 11 years (2005-2015) the validation period. Figure 2 shows the spatial distribution of the collected data used in this work. A brief description of these data follows, as well as the setup of the SWAT project.

Meteorological Data
Due to the low spatial-temporal resolution of observed data pertaining to temperature, solar radiation, relative humidity, and wind speed, this work uses gridded daily meteorological data obtained from the National Center for Environmental Prediction-Climate Forecast System Reanalysis (CFSR), with global atmosphere spatial resolution of around 38 km. The data for total daily precipitation was provided by the Brazilian National Water Agency (ANA), which made available a collection of data from 149 institutions. As shown in Figure 2a, the study area has a good spatial density of stations, with 2494 rain gauges within the basin, with the majority located in the eastern side of the UPRB. The rain gauges have different data availability during the simulation period, that may range from only a few years of data, up to the total period of simulation with no missing data. About half of the rain gauges (47%) contain less than 20% missing data.
The precipitation data was thoroughly controlled before use. First, quality checks, such as double records, typos, and the location of stations were evaluated. Routines to assess such inconsistencies were developed in R programming language. After that, the inverse distance weighted (IDW) method [27] was used to interpolate the daily precipitation records over the basin to a resolution of 0.1 degree. IDW has been used in several studies for interpolation of precipitation over hydrological basins and provided satisfactory results [28]. Water 2019, 11, x FOR PEER REVIEW 5 of 20

Topography
The topographic features ( Figure 2b) were characterized according to a Digital Elevation Model (DEM) map at 30-m resolution obtained from the Shuttle Radar Topography Mission, available from http://srtm.csi.cgiar.org/srtmdata/. Based on this model the digital river network, as well as the sub-basins, were generated.

Soil Data
The soil map was elaborated from the information provided by the Brazilian Agriculture Research Corporation (EMBRAPA, 2011) at a scale of 1:5,000,000. For the Paraguayan portion of the basin, the Harmonized World Soil Database (HWSD, 2011) with spatial resolution of 1 km was used. The initial classification considered 25 classes of soil types. In this study, the characteristics of oligotropic, mesotropic, eutropic, and dystropic soils were grouped in a single class, resulting in 15 classes. The Dark-Red latosols and argisols are predominant in the basin (Figure 2c), representing 43.9% and 20% of the area. The properties of each soil class were collected from a diverse set of documents that used the SWAT model in Brazilian basins [29][30][31]. The list of soil parameter values adopted in the simulation can be found in the Supplementary Materials (Table S1).

Land Use and Land Cover Data
Land use and land cover data were obtained from the Rudke [32] classification. The original classification of 10 different categories was reclassified into six dominant classes according to the SWAT land use classification ( Figure 2d). As a result, the Agriculture Land-Generic (AGRL) and Pasture (PAST) are the main classes and comprise 46.1%, and 25.6% of the total area, respectively. They cover mainly the western portion of the basin. The next two major classes are Forest-Evergreen (FRSE) and Range-Grasses (RNGE), encompassing 20.2%, and 5% of the basin, respectively. The remaining area of the basin is covered by Water (WATR), which covers 2%, and Residential Med/Low Density (URML), 1.1%. Most of the urban areas are concentrated in the headwaters of the main tributaries of the basin, such as in the upper Tietê and Iguaçu rivers.

Model Set Up
Based on the previously described data, the UPRB was discretized into 5,187 sub-watersheds, using a threshold drainage area of 100 km 2 , with an average size of about 173 km 2 (see Figure 2e). For most applications, the default threshold values used to define HRU's are 20%, 10%, and 20%, for land use, soil type, and slope, respectively [33]. However, in order to allow the assessment of land use changes in future research, further details are needed. Hence, the resulting sub-watershed was defined by the combinations of land use, soil types, and slope, using a threshold of 5%, 10%, and 20%, respectively. As a result, 44,635 HRU's were generated. In addition, five categories of slope were defined as this is the maximum number of categories possible. They are flat (0-3%), smooth rolling (3-8%), wavy (8-20%), strong wavy (20-40%), and hilly (>45%), according to EMBRAPA [34] classification. The potential evapotranspiration (PET) was calculated using the Penman-Monteith method [35] and the surface runoff within the model with the Soil Conversation Service's Curve Number method [36].

River Discharge Data
Monthly river discharge data were organized based on calibration period  and validation period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The data comprise both natural streamflow data, provided by theNational Water Agency (ANA), and naturalized discharges, obtained from the National Electrical System Operator (ONS). Only discharge series with at least 32 years of daily records and less than 20% of missing data were selected. For the western side of the basin, however, a threshold of 40% of missing data was used due to the low quality of the data available.

SUFI-2 and Parameters Calibration
Calibration, sensitivity, and uncertainty analysis were performed by the Sequential Uncertainty Fitting (SUFI-2) algorithm proposed by Abbaspour et al. [37], using SWAT-CUP version 5.1.6.2 [38]. Moreover, to optimize the model execution, the parallel processing module [39] was used. SUFI-2 was developed by considering the uncertainties of parameter ranges, which are sampled through Latin hypercube sampling. The main aim of this algorithm is to estimate the most observed variables within the 95 PPU band, which is quantified at the 2.5% and 97.5% of the cumulative distribution. SUFI-2 considers two indices to evaluate the performance of the calibration: The p-factor, calculated through the percentage of observed variable bracketed by the 95 PPU, which varies (between 0 for useless simulation and 1 for perfect simulation); and the r-factor, calculated through the ratio of the average width of the 95 PPU band (Prediction Uncertainty) and the standard deviation of the observed variable.
SWAT contains a large number of parameters that describe the processes in the soil-plantatmosphere interface. To calibrate the discharge series, a list of 20 parameters related to stream flow was selected as shown in Table 1. The choice of parameters, as well as their ranges, was based on previous research [6,[40][41][42]. Parameters that govern the soil, SCS runoff curve number, soil available water storage capacity, saturated hydraulic conductivity, and soil evaporation compensation factor were used. The selected parameters that govern groundwater were: Threshold depth of water in the shallow aquifer for return flow, groundwater delay, threshold depth of water in the shallow aquifer for "revap", deep aquifer percolation fraction, groundwater "revap" coefficient, base flow alpha factor, and base flow alpha factor for bank storage. For the channel, effective hydraulic conductivity in channel and Manning´s value for the main channel were chosen. For the parameters governing

SUFI-2 and Parameters Calibration
Calibration, sensitivity, and uncertainty analysis were performed by the Sequential Uncertainty Fitting (SUFI-2) algorithm proposed by Abbaspour et al. [37], using SWAT-CUP version 5.1.6.2 [38]. Moreover, to optimize the model execution, the parallel processing module [39] was used. SUFI-2 was developed by considering the uncertainties of parameter ranges, which are sampled through Latin hypercube sampling. The main aim of this algorithm is to estimate the most observed variables within the 95 PPU band, which is quantified at the 2.5% and 97.5% of the cumulative distribution. SUFI-2 considers two indices to evaluate the performance of the calibration: The p-factor, calculated through the percentage of observed variable bracketed by the 95 PPU, which varies (between 0 for useless simulation and 1 for perfect simulation); and the r-factor, calculated through the ratio of the average width of the 95 PPU band (Prediction Uncertainty) and the standard deviation of the observed variable.
SWAT contains a large number of parameters that describe the processes in the soil-plant-atmosphere interface. To calibrate the discharge series, a list of 20 parameters related to stream flow was selected as shown in Table 1. The choice of parameters, as well as their ranges, was based on previous research [6,[40][41][42]. Parameters that govern the soil, SCS runoff curve number, soil available water storage capacity, saturated hydraulic conductivity, and soil evaporation compensation factor were used. The selected parameters that govern groundwater were: Threshold depth of water in the shallow aquifer for return flow, groundwater delay, threshold depth of water in the shallow aquifer for "revap", deep aquifer percolation fraction, groundwater "revap" coefficient, base flow alpha factor, and base flow alpha factor for bank storage. For the channel, effective hydraulic conductivity in channel and Manning´s value for the main channel were chosen. For the parameters governing land use and land cover factor, plant uptake compensation factor and maximum canopy storage were selected.
Finally, for the parameters governing the sub-basin, surface runoff lag time, average slope length, lateral flow travel time, and average slope steepness were used. * "r_" refers to a relative change in the parameters where the current values is multiplied by 1 plus a factor from the given parameter range.

Objective Function
To assess the performance of the model, it is recommended that the simulation should be evaluated by several statistical indices [23]. Five indices were chosen so that they, together, can provide a general overview of the quality of the simulations ( Table 2). The percent bias (PBIAS) [43] provides a measure of how consistently simulated values are higher or lower than observed ones. The coefficient of determination (R 2 ) provides the proportion of the variance of the original data that is explained by the simulated values. Nash-Sutcliffe efficiency (NSE) [44] measures whether an observed value is better estimated by the model result or by the average of observed values. R 2 generally enhances the fitting of the model to lower values, while NSE tends to emphasize the fitting of high values. The Kling-Gupta efficiency (KGE) [45] intends to be a more general index that compares the variability of the observed and estimated values by including information about the correlation between them and their standard deviations, as well as any bias present, which is expressed by the relation between the mean values. Finally, the ratio between the standard deviation of the observations and the root mean square error (RSR) [46] measures how large the standard error is compared to the variability of the original data. Table 2 presents the equations of the objective functions used in this work and their model performance rating based on the threshold suggested by Moriasi et al. [46] and Thiemig et al. [47] for monthly discharge. These works classified the simulation into four performance ratings: Uunsatisfactory, satisfactory, good, and very good.
Ratio of standard deviation of observations to root mean square error (RSR)

Modelling Protocol
The criteria and the procedures used for the calibration and validation processes are summarized as follows: I.
In order to run the simulation with parallel processing, due to memory limitations as a result of the project size, the basin area was divided into 9 watersheds for calibration and the fitted values in each sub-basin were used for the initial project. II.
The geographic position of each outlet was verified. According to previous modelling studies [6,48], one of the main calibration problems is the incorrect position of the outlets. III.
A multi-objective calibration, which consists of simultaneous multi-site calibration from upstream to downstream outlets, was performed. This technique was recommended by Leta et al. [49] for a heterogeneous basin and presented better results compared to other methods such as single-site calibration (SSC). IV.
The discharge outlets which performed satisfactory or better in all objective functions that are presented in Table 2 were not considered in the calibration process. V.
The initial parameter ranges followed the calibration protocol presented by Abbaspour et al. [48] for large-scale basins. For example, if the simulation presented base flow too low (high), the GWQMN, GW_REVAP, and REVAMPM parameters should increase (decrease). Therefore, before each calibration, the temporal evolution of the discharge simulation was evaluated as to whether it underestimated or overestimated the observation. VI.
SUFI-2 provides several objective functions for calibration. The objective function selected in the calibration process was NSE. This index has been used in several studies and provided satisfactory results [50]. VII.
Once the sub-project was built for the sub-basin, and the ranges of parameters were defined, the model simulations were run between 150 and 500 times, with a maximum of 3 iterations. The numbers of simulations, as well as of iterations, were based on the size of the sub-project and performance of the initial simulation.

Sub-Basins Selected for Calibration
As stated in the modelling protocol, the criterion for selecting outlets for calibration was that they had an unsatisfactory rating in at least one of the objective functions presented in Table 2 over the simulation period . Figure 4 illustrates the sub-basins used for calibration (top figure), and examples of the temporal evolution of observed and simulated monthly discharge, as well as the values of performance indicators. Considering the 78 outlets selected, 23 performed satisfactory or better in all objective functions following the classification suggested by Moriasi et al. [46] and Thiemig et al. [47]. Hence, these sub-basins outlets were not used for the calibration process. It is clear that most of the outlets that performed well are located in the southern parts of the basin, especially in the Iguaçu sub-basin (VI) and adjacent areas of the Paraná (IV) and Paranapanema (V) sub-basins. This goodness-of-fit between measured and simulated discharge is mainly due to a large number of precipitation stations that are located over the sub-basins, with a low percentage of missing data (see Figure 2). For instance, in the streamflow of the Upper Iguaçu River (Figure 4b), the model has a good representation of the average, minimum and maximum discharge values. Regarding the statistical indices, SWAT has provided more than satisfactory results with 7.8 (very good), 0.77 (very good), 0.74 (good), 0.86 (good), and 0.51 (good) for PBIAS, R 2 , NSE, KGE, and RSR, respectively. The remaining figures of the temporal evolution of the outlets that yielded good performance are available in Figure S1 in the Supplementary Materials.

Calibration and Validation Performance
Figures 5a to 5e show the spatial distribution of the values of the objective functions used to evaluate the goodness-of-fit of measured discharge data estimated by SWAT. The performance of the monthly simulations for the calibration  and validation (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) period ranged from very good to unsatisfactory. It is clear that after the calibration process, the model has a good representation of monthly discharge values for most of the outlets of the UPRB that are located mainly in the Grande (II), Tietê (III), Paranapanema (V), and Iguaçu (VI) sub-basins. On the other hand, Paraná (IV) and Paranaiba (I) were the sub-watersheds that had the highest number of outlets with unsatisfactory simulations. This can be attributed to the low density of rain gauges mainly on the Ivinheima and Sucuriú river basins located on the Paraná sub-basin (IV), on the western side of the basin.
The indices R 2 ( Figure 5a) and PBIAS (Figure 5b) present the best hydrological performance for all sub-basins, with 92% and 86% of the outlets showing satisfactory or better performances. For R 2 , 61 (78%) of the outlets performed better than satisfactory with values of up to 0.91 over the Paraná and Sapucaí rivers. Similarly, the PBIAS index gave more than half of the outlets (63%) a better than satisfactory rating.
The rating of the KGE index (Figure 5d), which is based on the equal weighting of three different components (correlation, bias, and variability), shows that 88% of the outlets performed better or equally satisfactory. Only 9 outlets produced unsatisfactory simulations. The maximum value obtained for KGE was for the Grande river with 0.87.
Finally, the NSE ( Figure 5c) and RSR (Figure 5e) indices were those with the highest number of outlets with unsatisfactory simulation. However, the percentage of satisfactory stations was still high. Considering NSE > 0.5 or RSR < 0.7 for a satisfactory simulation, the model reached this criterion in 76% and 74% of the outlets, respectively. One of the reasons that explain why these indices performed slightly below the others is the low quality of the simulations of the base flow. This limitation is underlined by previous studies that evaluated the hydrological routines of the SWAT model [51]. SWAT simulates two types of aquifers: shallow (unconfined) aquifers, which contribute to return flow to streams within the catchment, and deep (confined) aquifers, which are responsible for the flow outside the basin (amount of water used, for example, for irrigation and water supply) and are considered water sinks in the system [26]. Once the model calculates the groundwater, studies that present difficulties in representing transfers associated with these types of water may present an  The indices R 2 ( Figure 5a) and PBIAS (Figure 5b) present the best hydrological performance for all sub-basins, with 92% and 86% of the outlets showing satisfactory or better performances. For R 2 , 61 (78%) of the outlets performed better than satisfactory with values of up to 0.91 over the Paraná and Sapucaí rivers. Similarly, the PBIAS index gave more than half of the outlets (63%) a better than satisfactory rating.

Calibration and Validation Performance
The rating of the KGE index (Figure 5d), which is based on the equal weighting of three different components (correlation, bias, and variability), shows that 88% of the outlets performed better or equally satisfactory. Only 9 outlets produced unsatisfactory simulations. The maximum value obtained for KGE was for the Grande river with 0.87.
Finally, the NSE ( Figure 5c) and RSR (Figure 5e) indices were those with the highest number of outlets with unsatisfactory simulation. However, the percentage of satisfactory stations was still high. Considering NSE > 0.5 or RSR < 0.7 for a satisfactory simulation, the model reached this criterion in 76% and 74% of the outlets, respectively. One of the reasons that explain why these indices performed slightly below the others is the low quality of the simulations of the base flow. This limitation is underlined by previous studies that evaluated the hydrological routines of the SWAT model [51]. SWAT simulates two types of aquifers: shallow (unconfined) aquifers, which contribute to return flow to streams within the catchment, and deep (confined) aquifers, which are responsible for the flow outside the basin (amount of water used, for example, for irrigation and water supply) and are considered water sinks in the system [26]. Once the model calculates the groundwater, studies that present difficulties in representing transfers associated with these types of water may present an unsatisfactory performance for the base flow prediction with the SWAT model. For instance, Srivastava et al. [52] found a NSE value of −0.16 in the predictions of monthly base flows. Similarly to the current study, Wu and Johnston [53] simulating long-term periods found it difficult to simulate dry seasons with the model. In this case study, the SWAT model performed better in simulating wet seasons than dry seasons. unsatisfactory performance for the base flow prediction with the SWAT model. For instance, Srivastava et al. [52] found a NSE value of −0.16 in the predictions of monthly base flows. Similarly to the current study, Wu and Johnston [53] simulating long-term periods found it difficult to simulate dry seasons with the model. In this case study, the SWAT model performed better in simulating wet seasons than dry seasons.   Figure 6 shows the comparison between the observed data and simulated values for the temporal evolution of the monthly discharge in the calibration and validation  period. The plots show the final outlets of the main rivers of the UPRB. Even though the model did not have a good estimate of the discharge at some outlets in the basin, these did not have a significant effect on the final outlet of the main rivers of the basin, due to their contribution area. This could be explained by the difference among the magnitudes of discharges. For instance, a closer examination of the long-term monthly mean discharge at the final outlets of the rivers shows that the Paranaíba river has a discharge of 2465 m 3 s −1 , while the Da Prata river, one of its tributaries has an average discharge about 71 m 3 s −1 , which represents 3% of Paranaiba river. The fact that the simulation for the Da Prata river performed an unsatisfactory simulation in R 2 , NSE, and RSR indices did not impact the quality of the performance of the Paranaíba. Similar cases occur in other major rivers of the UPRB. Figure S2, in the Supplementary Materials, shows the remaining graphs of the temporal evolution of the discharge on outlets after the calibration process.
Water 2019, 11, x FOR PEER REVIEW 14 of 20 Figure 6 shows the comparison between the observed data and simulated values for the temporal evolution of the monthly discharge in the calibration and validation  period. The plots show the final outlets of the main rivers of the UPRB. Even though the model did not have a good estimate of the discharge at some outlets in the basin, these did not have a significant effect on the final outlet of the main rivers of the basin, due to their contribution area. This could be explained by the difference among the magnitudes of discharges. For instance, a closer examination of the long-term monthly mean discharge at the final outlets of the rivers shows that the Paranaíba river has a discharge of 2,465 m 3 s −1 , while the Da Prata river, one of its tributaries has an average discharge about 71 m 3 s −1 , which represents 3% of Paranaiba river. The fact that the simulation for the Da Prata river performed an unsatisfactory simulation in R 2 , NSE, and RSR indices did not impact the quality of the performance of the Paranaíba. Similar cases occur in other major rivers of the UPRB. Figure S2, in the Supplementary material, shows the remaining graphs of the temporal evolution of the discharge on outlets after the calibration process.   Table 3 shows the objective function values from the final outlets of the main rivers for the calibration  and validation (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) periods. PBIAS ranged from satisfactory to very good simulation both for calibration (mean = −7.86) and validation (mean = −15.5) for the five rivers. High values of R 2 , greater than 0.80, were found in both calibration and validation results, indicating a very good correlation between the monthly observed and simulated discharges. In the calibration period, the NSE and KGE ranged from 0.56 (satisfactory) to 0.73 (good), and from 0.55 (satisfactory), to 0.77 (good), respectively. In the validation period, the NSE and KGE ranged from to 0.51 (satisfactory) to 0.73 (good), and from 0.55 (satisfactory), to 0.67 (satisfactory), respectively. Regarding the RSR index, values between 0.52 (good) and 0.66 (satisfactory) were found during the calibration process. For the validation process, only the Paraná river represented an unsatisfactory simulation with RSR value of 0.71. The remaining objective functions values of the outlets over the UPRB can be found in the Supplementary Materials (Table S2). As a whole, the calibration and validation of the outlets of UPRB provided promising results as indicated by acceptable values of statistical indices. The performance is better or comparable to other SWAT applications over Brazilian watersheds. For instance, Creech et al. [54] reported NSE ranging from 0.42 to 0.75 and from 0.42 to 0.77 for monthly discharge calibration and validation periods of the São Francisco River, the largest basin in the northeast of Brazil. On the other hand, considering small basins, Rocha et al. [55], modelling São Bartolomeu Stream Watershed, showed values of NSE and R 2 indices between −1.19 and 0.91, and 0.22, and 0.96, respectively. In addition, the results presented here agree with the range found in previous works where SWAT was calibrated for large basins worldwide. For example, Pagliero et al. [56] estimated the monthly flow for representative regions of the Danube basin found NSE ranging from 0.22 to 0.75 and R 2 ranging from 0.68 to 0.88. Another study performed by Easton et al. [57] for the Upper Blue Nile Basin showed values of R 2 ranging from 0.73 to 0.92 and NSE from 0.53 to 0.92.
One of the strengths of the current work is that the simulation was performed at a high spatial resolution, with the basin being divided into 5187 sub-basins and further into 44,635 HRUs. In addition, the project was built for a long-term simulation over 37 years . These spatial and temporal resolutions were not found in previous studies of large-scale SWAT applications. For instance, Jha et al. [58] simulated the streamflow of the Upper Mississippi River, which has an area around 447,500 km 2 , and discretized the basin into 119 sub-basins. These represent around 3760 km 2 of the average sub-watershed area, compared 179 km 2 for the current study basin. Pagliero et al. [56] defined 4663 HRUs (10% of our HRUs basin) over the Danube Basin, which has a drainage area of about 803,000 km 2 .

Summary and Conclusions
In this work, the Upper Paraná River Basin was built with the highest possible spatial discretization using the SWAT model for a long-term period between 1979 and 2015. The model was calibrated and validated using the SUFI-2 method for a large number of outlets. In addition, the evaluation of the performance of the model was carried out using several objective functions. The following conclusions can be drawn: The methodology used in this work regarding data preparation, model setup, and strategies for calibration and validation, as well as evaluation can be used for other large scale basins, especially in South America. II.
Due to the high spatial resolution and the good quality of most datasets collected in both meteorological, and physical variables, 23 outlets over the basin performed satisfactory or better in all the objective functions evaluated without the calibration process. Most of these outlets were found in the Iguaçu sub-basin (VI). III.
After the calibration process, most of the outlets analyzed (≥74%) presented better or equally satisfactory in all objective functions, mainly in the southern basin, which is the region with the highest density of stations. IV.
Although there are outlets with some errors in the simulated discharge, most of the evaluated outlets in the basin are in agreement with the observation especially those located at the final outlet of the main rivers of UPRB, which have the most significant contribution for the final discharge of the basin project.
The results provided in this work could be used for evaluating the potential impacts of land use and land cover as well as climate shift scenarios, which is the forthcoming study of the authors.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/5/882/s1, Figure S1: Temporal evolution of the monthly discharge and their statistical indices values without calibration process, Figure S2: Temporal evolution of the discharge rivers and their statistical indices values after the calibration process, Table S1: Soil parameters values used in the model, Table S2