Impacts of land use and cover change on Paraíba do Sul whatershed streamflow using the SWAT model

The main purpose of this work was to analyze the impact of land use and cover change in the Paraíba do Sul basin using SWAT model. Two land use and cover scenarios were analyzed, one referring to the year 1986 and other to the year 2015. The results showed that the SWAT model adjusted very well to the basin after calibration and validation. In addition, it was observed that in the scenario of 2015 the long-term flow rate was lower in a large part of the study area, especially the sub-basins in the stretch of the Paraibuna river. It was noticed that the changes in the behavior of the permanence curve were more representative in the Paraibuna river and in the Preto river. In the first one, due to the decrease of the forest area, the drought conditions and minimum flows were intensified. In the second, water availability improved in the 2015 scenario. Finally, the results presented here along with the calibrated SWAT model can help the water resource managers in the Paraíba do Sul basin and future works that aim to evaluate the impacts of climate change on the streamflow.


INTRODUCTION
Since the enactment of the 9.433/97 Law in Brazil, the management of water resources, as long as related to decisionmaking, is based on the limits of basins, which are natural areas of great economic, social and environmental importance. And on these such important regions a huge amount of anthropic activities are held such as industries, agriculture, deforestation, wildfires, urbanization among others.
From this perspective, some authors (Deng et al., 2015;Fan & Shibata, 2015;Neupane & Kumar, 2015;Sajikumar & Remya, 2015;Blainski et al., 2017;Gashaw et al., 2018;Zhang et al., 2018;Pereira et al., 2016;Silva et al., 2018) highlight that the land use and cover changes, mainly due to anthropic actions, on a basin have an effect on variables inherent to the water balance, such as the soil infiltration rate, the evapotranspiration, surface runoff and the water bodies recharge. Putting them all together, it's of major importance for the management of water resources to comprehend the impacts of the change of land use and cover on the basin runoff characteristics. Nowadays, the most indicates and used tools, for this purpose, are hydrological models.
Hydrological modelling is a powerful tool that attempts to simulate the hydrological processes of a basin. An integral part of any hydrological model is the rainfall-runoff process, which aims to estimate the runoff variables from real meteorological data, such as precipitation, temperature, solar radiation, wind speed and relative humidity. The hydrological models can be classified in three categories, physics-based models, conceptual models and black box models (Sajikumar & Remya, 2015). This classification is based on the actual rainfall-runoff process physics representation; models physically based have a more accurate representation of the hydric process (demanding specific entering information), easing the quantification the impact caused by the input data. Whilst the two others do not present this advantage (Sajikumar & Remya, 2015).
Currently, within the huge number of hydrological models, the Soil and Water Assessment Tool (SWAT) is one of the most used in studies that research the impact of the land use and cover different scenarios under flow regimes. SWAT is a semi-distributed model, continuous in time and physically-based, that was developed by the Agricultural Research Service and Texas A&M University having one of its applications the simulation of the impact when changing the use and cover of the land in basins using different scales (Arnold et al., 2012).
In a national scope, more than a hundred studies were registered using the SWAT model, aiming to assist the management of water resources (Bressiani et al. 2015). In Paraíba do Sul basin, an outstanding work was done by Pereira et al. (2016), which calibrated and validated the SWAT in the Pomba River watershed, in different points and using different validation tests. The authors had satisfactory results and concluded that the model may be used with good estimations of the region's water balance. Blainski et al. (2017) applied the SWAT model seeking to evaluate the flow alterations due to the changes on land use and cover observed from 1957 until 2012, in Camboriú basin in the Santa Catarina State, southern region of Brazil. They noticed, from the daily flow permanence curve, that the flowing between Q 60 (flow overcome in 60% of the time) and the Q 90 (flow overcome in 90% of the time and which represents the drought condition) was lower in the land use and cover 2012's scenario in comparison to the 1957's one. Silva et al. (2018) used the SWAT model to analyze the impact of different land use and cover scenarios on the water regime of the sub-basin of the lower-middle course of São Francisco River, in the northeast region of Brazil. The analyzed scenarios were: S1 -conversion of pasture areas into caatinga (native vegetation); S2 -conversion of pasture areas for agriculture; S3 -conversion of pasture areas into exposed soil. Silva et al. (2018) noticed that the flow and the debris transport augmented substantially in the S3 scenario.
This study seeks to calibrate and validate the SWAT model in the Paraíba do Sul basin, and to evaluate the variations occurred in the streamflow between the land use and cover, considering the land use and coverage scenarios of 1986 and 2015, aiming to assist the planning and management processes of this region.

Study area
The study area was the basin of Paraíba do Sul River, located in the southeastern part of Brazil between 20°26'S and 23°38'S latitudes and the longitudes of 41°00'W and 46°25'W, it has a drainage basin of 62,074 km 2 . There is a huge population density, approximately 17,634,301 inhabitants within its limits and encompasses, either partially or fully, 39 municipalities in the state of São Paulo, 57 in Rio de Janeiro and 88 in the state of Minas Gerais, adding up a total of 184 municipalities (Associação Pró-gestão das Águas da Bacia Hidrográfica do Rio Paraíba do Sul, 2012; Instituto Estadual do Ambiente, 2014).
Since the gauging station most downstream of the Paraíba do Sul River, with a good time range and few fails was chosen the 58974000 station (Campos -Municipal Bridge), the limits used on this study was its contribution area of 55,620km 2 . The limits of the watershed area, as well as the Paraíba do Sul basin's can be seen in the Figure 1.
The SWAT was developed in the early 90's and it is a physically-based model, semi-distributed (space-time variation), continuous in time and capable of simulating long periods of time. In the model, the basin is divided in multiple sub-basins which are now divided in Hydrological Response Units (HRUs), those are homogeneous areas which exhibit the same types of soil, use and cover of land and declivity (Arnold et al., 2012). Categorizing sub-basins in HRUs increases accuracy and provides a more refined physical description.
SWAT performs the simulation of the hydrologic cycle in two phases: a land one and the other is cycle variables propagation in a frame -routing phase. The ground phase consists of describing the water movement through the soil to the main channel, while the second one to estimates the way water is transported through the channels (streamflow) (Brighenti et al., 2016).
The land phase of the hydrological cycle is estimated by SWAT based on the water balance in accordance with Equation 1 (Arnold et al., 1998).
In which SWt is the final amount of water in the soil (mm), SW0 is the initial amount of water in the soil (mm), t is time in days, Rday is the total daily precipitation (mm), Qsur is the total daily surface runoff (mm), Ea is the total daily evaporation (mm), Wseep is the amount of water that seeped and got stored throughout the whole day (mm) and, finally, Qgw is the daily return flow (mm).
In this paper, the surface runoff Qsur was obtained using the method of the curve number (CN) from Soil Conservation Service (SCS), which is one of the most used methods (Ghoraba, 2015). The curve number method predicts the surface runoff given a rainfall event using a retention parameter, in millimeters, which varies according to the type of land use and cover the soil properties and hydrological conditions (Ghoraba, 2015;Gashaw et al., 2018). The CN method is obtained by the Equation 2.
In which Qsur is the surface runoff, mm, Rday is the total daily precipitation, mm e S is the retention parameter, mm. The retention parameter is given according to Equation 3.
. , 1 000 S 25 4 10 CN In which S is the drainable water volume in the soil of a saturated thickness area unity (mm/day), CN is the curve number.
Once the surface runoff is calculated by the curve number method, the routing phase starts. Thenceforth, the surface runoff that ends up in the main river of a basin is predicted by the calculation of the surface runoff from each sub-basin on its upstream.
In SWAT, there are three methods to calculate the potential evapotranspiration (PET), those being he Penman-Monteith Method, the Priestley-Taylor one and the method of Hargreaves (Abiodun et al., 2018). This paper uses the first one.
Once the qualitative description of how the model works is presented, its application is presented hereinafter.

Input data and model application
In this work, the ArcSWAT in ArcGis 10.1 was used. The applied procedures, from the first step until the hydrologic simulation can be seen in Figure 2.
The first step of the model consists in the basin's delineation, of the streamflow and the obtainment of morphometric parameters of the basin to be studied. Thenceforth, it was made necessary the use of a digital elevation model (DEM) that covers the area to be studied. In this paper, it was chosen to use the product HydroSHEDS (Hydrological data and maps based on Shuttle Elevation Derivatives at multiple Scales) whose spatial resolution reaches 90 meters and was produced by the conservation science program of the World Wildlife Fund (WWF) jointly with the United States Geological Survey (USGS). This DEM is a product obtained by the correction of Shuttle Radar Topography Mission (SRTM) data, using filter techniques and drainage conformity (Lehner et al., 2006).
In the step of delineation of the watershed, the contribution area of the station 58974000 (Campos -Ponte Municipal), was divided in 94 sub-watersheds. Concluding the delineation step, the procedure of creation of HRUs commenced ( Figure 2). For that, two land use and cover maps were used, from the years  1999). At last three classes of sloping were employed, from 0% up to 20%, 20% up to 40% and steeper than 40%. After defining the input data of this step, it was possible to create the HRUs: 389 hydrologic response units were created overall.
By the end of the creation of the HRUs, the next step is to create tables with the values of climate variables temperatures, precipitation, relative humidity, solar radiation and wind speed ( Figure 2). To do so, data of eleven INMET meteorological stations from 1985 until 2014 were used. The location of those meteorological station and three ANA's gauging station that were used to calibrate and validate the model and the limits of the main sub-basins of Paraíba do Sul River are in Figure 3. Before entering the meteorological data of the project in SWAT, the WGEN weather generator was applied (Figure 2) resulting in a series of statistics based on the data of the meteorological data, which will be added to SWAT so the model can rectify the flaws of the data.

5/13
Finally, having done all the steps of basin delineation, HRUs creation and importing weather data, the SWAT model is ready to simulate variables as: streamflow, evapotranspiration, sediments production and transportation, nutrients concentration etc. In this paper, monthly variables were simulated between the years of 2004 and 2014, without taking into account the two years of the model warm-up, so it would be calibrated and validated. The procedures done in the steps of calibration and validation are better detailed in the following section.

Calibration and validation
To calibrate a model means to adjust the output so that the simulated values approximate to the observed values. To do so, some adjustments are done in the most sensitive parameters of the model. Whereas validation aims to compare what has been simulated by the model calibrated using observed values, it is important to stress that the calibration and the validation periods must be different. That being done, the model is ready to be applied with accuracy (Arnold et al., 2012) In this paperwork, only the variable monthly streamflow was calibrated and validated. The period of calibration was from 2006 up to 2009, whereas the validation was from 2010 until 2014.
The software used to calibrate and validate was SWAT-CUP, in which the algorithm Sequential Uncertainty Fitting (SUFI-2) was chosen among the available ones in the system. This one was chosen since it presents good results with few iterations (Yang et al., 2008, Gashaw et al., 2018. SUFI-2 is a semi-automatic calibration and uncertainty analysis algorithm that takes into consideration all the uncertainty sources, including the uncertainty in the conduction variables. This one is expressed by the algorithm through the distribution of 95% of the probability (95PPU), which is calculated using the percentages of 2.5% and 97.5% of the cumulative distribution of an answer variable.
To evaluate the adjustments on the used parameters range, two statistical indicators, the p-factor and the r-factor. The first one varies from 0 up to 1 and indicates the amount of simulated data that is within the 95PPU range, values closer to the unity indicate a large quantity of simulated data in the reliability range. The second one indicates the length of the 95PPU range, in that case, the smaller the better. Abbaspour (2015) recommends a p-factor > 0.70 and a r-factor < 1.5.
The calibration and the validation were carried out in 3 points (Figure 3), in the gauging stations of the National Agency of Waters (ANA) coded as: 58520000 (Sobraji -MG), 58585000 (Manuel Duarte -MG) and 58974000 (Municipal Bridge of Campos -RJ). The parameters and their respective ranges that were used, for the monthly calibration, can be seen in Table 1.
On average, five iterations of 500 simulations were done for each station. In each simulation, the SUFI-2 tests a combination of values from the selected parameters. At the end of each simulation, the values of each parameters that produced the best result among the countless tested combinations are presented to the users, amongst various results. Before starting a new iteration, the calibration range of the value of each parameters was reduced (Abbaspour, 2015) which, besides presenting the best adjustment, presented also a high sensibility (value -p < 0.05). For validation, the data period between 2010 and 2014 was used.
Various statistical tools indicate the performance quality of the hydrologic model. In this paperwork, the coefficient of Nash-Sutcliffe (NS), Equation 4, the coefficient of residual mass (PBIAS), Equation 5, and the determination coefficient (R 2 ), Equation 6.   Addendum: The meaning of the letters which are previous to the parameter in the column Parameters: "a"; SWAT-CUP will add to the parameter a value within the range, "r" will multiply the parameter (1+ a value within the range) and, finally, "v"; SWAT-CUP will replace the value of the parameter by one within the calibration range.

7/13
In which, n is the total number of observations, Qobs is the observed flow (at an instant 'i'), Qsim is the simulated flow (at an instant 'i') e obs Q is the average flow observed in the analysis period.
NS varies from -∞ up to 1, in which the value 1 represents the perfect adjustment between the simulation and the observation and for NS < 0, it means that the average flow of the observed values is better than the simulated results (Brighenti et al., 2016). According to the classification of Motovilov et al. (1999) values of NS > 0.75 indicates that the result is very good; 0.36<NS<0.75 indicates a satisfactory result and NS<0.36 indicates that the result is unsatisfactory.
At last, the determination coefficient (R 2 ) indicates the adjustment between he simulated and the observed hydrograph, R 2 close to 1 indicates a very good adjustment.

Land use and cover scenarios
The images of land use and cover were developed from a mosaic of 9 images of the satellites Landsat 5 TM (1986) and 8 OLI (2015) from the database of the U.S. Geological Survey (USGS). These images were classified in 6 classes of use and cover: forest, urban area, exposed soil, water, pasture and agriculture. All the process of classification was executed by the software ENVI  5.2 ( Figure 2) using algorithms that classificates the maximum likelihood, one of the most used methods of supervisioned classification. To minimize the inherent errors of the classification process, Figure 5. Results from the calibration and validation step.

8/13
the software DINAMICA-EGO was used to correct and refine the image classification in areas that the error of the algorithm's classification was noticeable. The validation of the classification was done by the kappa's coefficient. The maps of land use and cover from 1986's and 2015's scenarios can be found in Figure 4.
The land use and cover applied in the basis design of SWAT, in other words, the one that was calibrated and validated was the one from 2015. The map of 1986 was applied after the proper validation and calibration so the impacts of land use and cover on the streamflow could be examined in this period.

Calibration and validation of the model
The results of the calibration step (2006 until 2009) and validation (2010 until 2015) of the monthly streamflow (m 3 /s) for each one of the analyzed gauging stations can be seen in Figure 5. It is noticeable that both in the calibration period and in the validation one, the SWAT model exhibit a great adjustment, according to the classification written by Almeida et al. (2018). In which the model is considered great if 0.75 < NS < 1.00, |PBIAS| < 10 and 0.75 < R 2 < 1.00.
The streamflow rates were underestimated in all the calibrated gauging stations for they presented positive PBIAS . Only in the validation period the stations 58974000 (Municipal Bridge of Campos -RJ) and 58520000 (Sobraji -MG) had their streamflow rate overestimated by the model.
Generally, the analyzed performance rates (NS, PBIAS, R 2 , p-factor e r-factor) obtained in the calibration and validation periods indicates a very good performance rate ( Figure 5). Noticing that, it is possible to say that the model can be used to evaluate the impacts of changes on land use and cover between the 1986's scenario and the 2015's scenario, being said that the only possible alteration is in the map of land use and cover, keeping the other variables, such as constant climate variables (model forcing) (Zhang et al., 2018).

Contribution of the land use and cover changes to the variation of streamflow
Contribution of the land use and cover changes to the in Table 2, the percentages of area from each presented class of land use and cover that can be seen in this paperwork of the study area. It is noticeable that the major class of land use and cover is pasture, which increased in 7,9% between 1986 and 2015.
The second most expressive class is forest, which decreased 3,7%. In this context, Tucci (2001) explains that in areas of either few or none green cover, the processes of interception and absorption of water by the roots of the plants gets damaged and starts to occur in a less effective way. This state of the soil eases the surface runoff and damages the process of water infiltration in the soil variation of streamflow.
To evaluate the impact of land use and cover changes in monthly streamflow, a simulation of the period from 1987 until 2014, ignoring the two first years to warm-up the model. In the map of Figure 6, there is the variation of average flow in long period (QMLP), which is the average of the average annual streamflow for all the data simulated by SWAT for all the 94 sub-watersheds. It is noticeable that most of the times the QMLP decreased between the scenarios of 1986 and 2015, highlighting some sub-watersheds that are part of the watersheds of the Paraibuna Rivers (in the dark red in Figure 6), in which the QMLP decreased more than 10% of its former value due to alterations in land use and cover from 1986 until 2015. However, in the Preto River basin, had an increase between 0% and 2,5% of the QMLP.
Despite the long period flow, there was also a comparison between the minor (Qmin) and the major (Qmáx) monthly streamflow (m 3 /s) which occurred during the period from 1987 and 2015 for the two scenarios of use and cover. Taking into consideration Qmax, it can be seen in the Figure 7 that, generally, the streamflow peak increased in the sub-basins located in the region of middle Paraiba do Sul, and decreased more than 15% in great part of the sub-watersheds of Paraibuna River and Dois Rios River. In the sub-watershed of Piraí River the maximum streamflow diminished more than 20%.
In Figure 8 there are three boxplot graphics that represent the distribution of minimum annual streamflow for each of the three studied stations. In Manuel Duarte station (Figure 8a), which is in the Preto River, the alteration of minimum annual streamflow between land use and cover scenarios of 1986 and 2015 were almost unnoticeable being slightly greater in 2015. However, in the Sobraji Station (Figure 8b), which is in the Paraibuna River upstream of the municipality of Juiz de Fora, the minimum annual streamflow had a huge decrease between the 1986's scenario and the 2015's scenario. Furthermore, the simulation showed that the distribution of the minimum annual streamflow in the station of the municipal bridge of Campos -RJ (Figure 8c) which is in the Paraiba do Sul River, was lower in the 2015's scenario.
In spite of the mentioned approaches, it was also researched the evaluation of impacts in the change of land use and cover between the 1986's and 2015's scenarios through the permanence curve, which indicates the frequency of a streamflow to be equal or greater in a water course, from the points which coincide with the 58520000 stations (Sobraji -MG) in the Paraibuna River ( Figure 9B), 58585000 (Manuel Duarte -MG), in the Preto River ( Figure 9A) and 58974000 (Municipal Bridge of Campos -RJ) in the Paraiba do Sul River ( Figure 9C). This curve was divided into four parts: Maximum streamflow (0% up to 10%), intermediary streamflow (10% until 60%), drought conditions (60% until 90%) and minimum streamflow (>90%).  One notes that the conditions of monthly streamflow in permanence curve in Figure 9B and Figure 9C were more positive, considering the fact that there was more water availability in the 1986's scenario in contrast with the 2015's one. In Figure 9B, it's noticeable that the drought conditions and the minimum streamflow suffered a slight reduction compared with the 1986's scenario. However, as Figure 9C indicates, the streamflow in the condition of drought and medium flow in 2015's scenario were inferior to the 1986's scenario. Nevertheless, the opposite picture is portrayed in the permanence curve of Manuel Duarte (58585000), in which the 2015's scenario of land use and cover reported a better water availability rather than the 1986's one, as shown in the permanence curve in Figure 9.A. Table 3 presents the values, in percentage, of the areas of each land use and cover in the regions upstream from the stations 58585000 (Manuel Duarte -MG) and 58520000 (Sobraji -MG). the forest area dwindled and the grazing land area expanded between 1986 and 2015. Still, the opposite effect can be seen in the 58585000 (Manuel Duarte -MG), expansion of the forest area and reduction of the grazing one. Andréassian (2004) and Tucci (2001) emphasize that areas in which the forest area increases, there is a drop in the flow peak and a rise in the recharge of watercourses due to the influence of this type of soil cover in opposition to grasslands or areas with no vegetation at all.  The increase of forest areas can be linked to the fact that a huge number of restrict use areas, such as the conservation unities (UC) and Private Natural Heritage Reserve in the watershed region of the Preto River in comparison with the other sub-basins of Paraiba do Sul River, as presented in Figure 10. It's remarkable, in the watershed of Preto River, the Area of Federal Environmental Protection of Serra da Mantiqueira (Sustainable Use), which was created by the Decree 91.304/85 and magnified in 1995 by the Law 9.097, part of Itatiaia National Park (Fully Protected), which was created after the Decree 1.713/97 and had its area expanded after the Decree 87.586/82; and last of all, the Area of Municipal Environmental Protection Boqueirão da Mira (sustainable use), created after the Law 929/01.

CONCLUSION
It was possible to evaluate the impacts on streamflows due to the changes of land use and cover. In general, the water availability in the Paraíba do Sul basin decreased between the 1986's and the 2015's scenarios According to the simulated values, the most impacted areas because of changes on land use and cover were the sub-watersheds upstream the station 58520000 (Sobraji-MG), which constitute the watershed of the Paraibuna River. This one displayed in all performed analyses, a diminishment on its water availability.
The analysis of the permanence curve from 58585000 (Manuel Duarte -MG), on the Preto River, depicted a rise on the minimum flows and henceforth, an improvement of the drought conditions in the land use and cover scenario of 2015, in comparison with the 1986's one. Furthermore, this fact is most likely to be related to the creation and amplification, under legal devices, of the Conservation Unities (federal and municipal) and of the Private Natural Heritage Reserve, in the years of 1986 and 2015.
Considering the hydrological behavior of most of the basin of Paraiba do Sul River, represented by its station located closer to the basin mouth -station 58974000, Municipal Bridge of Campos (RJ) -it was noticeable that, in general, the water availability of the study area decreased between the 1986's and the 2015's scenarios. Table 3. Percentage of the classes of land use and cover for the contribution areas of the stations 58585000 (Manuel Duarte -MG) and 58520000 (Sobraji MG).

Classes
Station 58520000