Effects of Deforestation on Spatio-Temporal Runoff Patterns in the Upper Teles Pires Watershed , Mato Grosso , Brazil

This study is an exploratory, spatio-temporal analysis of river discharge in the upper Teles Pires contribution area for the period between 1993 and 2006, and its relationship with deforestation dynamics on a watershed scale. After an evaluation of temporal variability of monthly rainfall and runoff time series, generalized linear models (GLM) have been fitted to examine the relationships between discharge with indicators of deforestation in the contribution area of six subwatersheds. Monthly discharges at the main outlet show a significant positive trend, whereas a slight increase of precipitation has no statistical significance. Besides precipitation, principal determinant of specific discharge, deforestation percentages and the annual relative rates of deforestation are significant predictors in the three smaller subbasins with areas up to about 14.000 km. In the most altered sub-basin under intensive agricultural use, the ratio of residence and hydrological response have positive correlations with deforestation indicators, evidencing the removal of arborous formations diminish evapotranspiration and favor discharge. With increasing contribution areas up to about 54,000 km, the interference of deforestation on runoff and hydrological parameters decrease and precipitation remains the unique significant predictor. Therefore state authorities should extend the monitoring and foster management in smaller watersheds, where negative impacts of deforestation on runoff were found to be more expressive than in the great, already monitored streams.


Introduction
Deforestation of savannahs and forests for crop and cattle farming is considered one of the major anthropogenic factors on the water balance, river discharge and water resource availability (Sahin & Hall, 1996;Bosch & Hewlett, 1982).Vegetation cover removal tends to reduce storage by interception, litter and in the root zone, changing groundwater recharge principally during strong precipitation events (Bronstert et al., 2002).Evapotranspiration is the main process causing changes in annual discharge after deforestation (Brown et al. 2005).Generally, the reduction of evapotranspiration triggers an increase in annual discharge.This effect, however, varies in different climates, under different physiographic characteristics of the watersheds and along the intra-and interannual cycles (Hibbert, 1967;Andreássian, 2004;Zhang and Schilling, 2006).Discharge increases after deforestation are stronger during the rainy season, favoured by a reduced infiltration due to a degradation of soil structure (Bronstert et al., 2002;Zuazo and Pleguezuelo, 2008).During dry periods runoff may decrease, induced by channel sedimentation and reduced retention by soil loss (Harr et al.,1982;Cheng, 1989;Cornish, 1993;Gustard and Wesselink, 1993).
The influence of deforestation on discharge was found to be non-stationary.Stronger short-term effects may be tampered after some years or decades (Bosch and Hewlett, 1982) when runoff can achieve new equilibria (Brown et al. 2005).
Hydrological alterations after Land Use and Cover Change (LUCC) are overlaid by climatic variability (Frans et al., 2013) and are scale-dependant (Blöschl et al., 2007).Discharge alterations due to LUCC are frequently local phenomena and are buffered in great basins (Woods, 2005;Defries & Eshleman, 2004), if compared with the effects caused by climatic variations.Whereas deforestation of small watersheds (<100 km 2 ) was found to systematically increase discharge, such effects cannot be generalized for the case of large basins (Bruijnzeel and Sampurno, 1990).In tropical sub-humid tropical watersheds, Costa et al. (2003) attributed increase in average and maximum discharge of the Tocantins river to the increase of agricultural land use, whereas Gyau-Boakye and Tumbulto (2006) observed a reduced discharge after substitution of forests by pasture and crops.
On the background of the scarcity of studies on LUCC on discharge in neotropical savannahs and adjacent transitional forests region (Oliveira et al., 2015) and contradicting results on the effects of deforestation on discharge (Krishnaswamy et al., 2012;Andreassian, 2004), we examine spatio-temporal variation of river runoff on the watershed scale in the upper Teles Pires river contribution area, which includes portions of one of the most productive crop farming region worldwide and major projects for the installation of new hydropower plants.

Study area
The upper Teles Pires watersheds drains central and northern parts of the state of Mato Grosso, between the headwaters in the Serra Azul mountain ranges and downstream its confluency with the Peixoto de Azevedo river between about de 10° e 14°50' lat.S. e 54° e 56°40' lon.W. (Figure 1).Zeilhofer, P., Alcantara, L. H., Fantim-Cruz, I. Time series of the outlets from six subwatersheds have been analysed, which have contribution areas between 5.491 km² and 52.407 km².Five stations are hereby gauging the Teles Pires and one the Rio Verde, its principal tributary.
Average elevations and slopes are varying between 462 mNN (Porto Roncador) and 300 mNN (Indeco) and 7,1° (Porto Roncador) and 4,1° (Cachoeirão), respectively (Table 1).The tropical semi-humid Aw climate has a major gradient of average annual precipitations between about 1500 mm in the south up to 2100 mm in the northwest.In the sub-watersheds, annual rainfall varied between 1,747 mm in the Rio Teles Pires headwater region to 1,829 mm in the entire watershed contribution area (17340000).
Specific discharge varied between 667 mm in the Rio Verde contribution area (17230000) and 798 mm at the second TP river gauge the higher part of the watershed yet (17210000).Highest specific discharge and interannual variations occur in the Teles Pires headwater regions (17200000,17210000), where terrain slopes are steeper.
Headwaters with elevations up to 900 mNN are located in residual mountain ranges of the Upper Paraguay and the Interplanaltic Depression of Paranatinga, where shallow Cambisols and Arenosols prevail.The middle parts of the watershed in elevations between 350 mNN and 550 mNN are formed by sedimentary formations of the Parecis highlands.Here, deep Ferralsols (Latosols) develop, sustaining intensive cash crop farming.The lower Teles Pires passes igneous and metamorphic rocks in a rolling relief with Acrisols/ Lixisols/Alisol of varying soil depths, mainly used for cattle farming.

Watershed parametrization
Sub-catchment contribution areas were delimited and physiographically parametrized by the 3 arc-seconds (90m-resolution) SRTM digital elevation model (DEM), version 4, obtained from USGS.For the delimitation procedure only, the original DEM was corrected by the TOPOGRID algorithm (Hutchinson, 1989), with imposition of a single line hydrographic network in a 1:100.000scale, obtained from the Mato Grosso State Environmental Agency (SEMA).Slopes and elevations were extracted from the non-modified DEM.Deforestation data for the period between 1992 and 2007 were obtained from SEMA.Time series gaps at the years of 1996, 1998 e 2000 were filled in by linear interpolation.

Hydrological data and pre-processing
Hydrological data sets of daily rainfall and discharge from January 1993 through December 2006 were obtained from the Hidroweb website maintained by National Agency of Water (ANA).Data gaps in both series were filled in by linear regression using data from the most highly correlated adjacent stations.From the 57 rainfall gauges in the watershed and adjacent areas, those 21 with less than 30% of data gaps were maintained.The six river gauges had gaps between 4 and 41 %, which were filled in by linear regression as well.The average coefficient of determination (R²) was 0.84 for gap filling in the daily discharge time series.For comparisons with precipitation and deforestation, monthly and annual discharges were determined.
Monthly areal precipitations by catchment were estimated by the Thiessen polygons method, Zeilhofer, P., Alcantara, L. H., Fantim-Cruz, I. which has been found to perform well if, compared with more complex interpolation procedures such as krigging, in slightly rolling relief (Burrough and McDonnell, 1998;Siska et al., 2005;Bertoni and Tucci, 2004).
For comparison with precipitation and calculus of the runoff ratio (RR), discharge data were transformed in specific runoff.
Evapotranspiration (ET) was then estimated by the difference between precipitation and discharge (Eq.1).
For comparison with annual deforestation, we calculated as well the Runoff Ratio (RR).RR is a dimensionless index, obtained by the division of the average discharge volume and average rainfall volume over the monitored period (Ochoa-Tocachi et al., 2016) and describes the runoff generation in a watershed (Eq.2).

Time series analysis
Hydrological time series typically show serial autocorrelation, which bias outcomes in trend tests (Kouraev et al., 2007).Thus, we first determined the empirical autocorrelation functions to choose an adequate trend test (Dunn, 2005).Accordingly monthly precipitation and runoff data were tested for monotonic trends using the seasonal Mann-Kendall test (MKs) (Hipel and Mcleod, 1994), whereby the series is divided in temporal subsets of same períods of the year (ex.months).The MK statistics (Sk) and its variance for every block k is given by (Eq.3): Where xi is the time at observation i, yi the value of the variable at observation i, m the number of periods (ex.months) and nk the number of observations during period k.Davison and Hinkley (1997) recommended block bootstrapping for significance testing in trend extraction, where samples of the original data set are randomly and repeatedly subsetted.The tau statistics of the MK test of each permutation are hereby compared with the tau of the original set.MK with bootstrapping has been applied successfully in hydrological time series analysis (Douglas et al., 2000;Burn and Hag Elnur, 2002).
As precipitation is the main driver in runoff generation, rainfall trends can mimic or exaggerate possible land use and cover changes effects on runoff.Therefore we analysed precipitation time series as well, using them as a proxy of climate induced variations on runoff (Hung and Zang, 2004).
Deforestation data were obtained from the Mato Grosso State Environmental Agency (SEMA), which elaborates data layers since 1992 by visual interpretation of Landsat (E)TM data.For comparison with specific runoff and RR, total deforested areas (DA) and yearly rates of deforestation (YRD) where determined for gauge stations contribution areas.
Generalized linear models (GLM) were applied to quantify relationships between precipitation, deforestation totals and yearly rates of deforestation with specific runoff and RR.We opted for GLMs, as they do not demand linearity between predictors and dependants, neither homogeneity of variance of the dependant variables (Nelder & Wedderburn, 1972).Allow tests of predictive power of ordinal-scaled data and the construction of nested models.Formal descriptions can be found in Dobson (1990).
In the six sampled watersheds, monthly runoff shows highest correlations with the precipitations of the previous month (Figure 2).Thus, annual discharges were compared with annualized precipitations with a one-month delay.Five of the six gauges are monitoring the Teles Pires river.Thus, their observations cannot be considered independent.Therefore, GLMs were generated separately for each contribution area.
Figure 2. Pearson correlation between monthly runoff (Qm) and rainfall without (Pm) and with time lags of one (Pm-1) and two (Pm-2) months for the six gauged contribution areas.
Goodness of fits of the GLMs was realized by the Likelihood Ratio Chi-square (LRC) and determined the probability (p<0.05) of the difference to the null model.Relative contribution of individual coefficients to the GLM was evaluated by the Wald Chi-square (χ 2 ) and its significance test (Dobson, 1990).All statistical analyses were realized, using the R package, version 3.2.1 (Cran).

Dynamics of deforestation, rainfall and discharge
Before the beginning of systematic monitoring by the State Environmental Agency (FEMA, actually SEMA) in 1992, already 13,776 km 2 or 26.3% of the basin had been deforested (Figure 3).
In 2007, 56.2% of the forest/woodland cover had been removed.Only the contribution areas of the gauge stations Roncador (17200000) and Teles Pires (17210000) in the upper basin had in 2007 more than 50% of natural vegetation preserved, yet, both having a stronger rolling relief and soils with low crop farming suitability (Tab.1).The Verde River watershed ( 17230000) is localized, almost entirely, in flat sandstone formations where Latosols (Oxisols) develop, which are intensively used for cash crop farming since the 1980ties.Here, a deforestation of 50% has been reached already in 1995.The Cachoeirão contribution area (17280000) had been deforested in 50% in 2000, whereas Tratex (17300000) and Indeco (17340000) achieved this rate in 2002.The mid and northern parts of the basin included igneous and metamorphic rocks with steeper slopes limiting deforestation pressures (Table 1).
YRD followed a single pattern throughout the contribution areas with the maxima in 1995 with rates higher than 6% for the watersheds 17200000, 17210000 and 17230000, 5.36% for 17280000, 4.79% for 17300000 and 4.15% for 17340000.A secondary peak had occurred between the years of 2002 and 2005, reflecting the economic grow and interest reductions in Brazil.In the southern part of the watershed, deforestation was proportional with the demand for crop farming, whereas the northern part was occupied by cattle farming, following in both parts an intense exploration of native woodlands and Since 2005, deforestation curves turned to be asymptotic.
Despite of the state of Mato Grosso has been already judged to control efficiently deforestation (Nepstad et al., 2014), we suppose that diminishing deforestation taxes occurred principally due to the reduction of non-protected areas with aptitude for agricultural production.Far parts of the basin in private areas are already deforested up or passing environmental legislation, causing a gradated reduction in deforestation (Soares et al., 2014).With its deforestation totals, the basin can be considered strongly altered, which has been in 2007 about 10% stronger deforested than averages in the entire Cerrado Biome (BRASIL, 2015).

Temporal variations of precipitation and runoff
Interpolated precipitations for the Indeco gauge (17430000) contribution area and the respective monthly runoff at the outlet have unimodal seasonal cycles (Fig. 4a, b), showing the typical hydrological regime of Tropical rainforest (Af) and Tropical Savannah (Aw) climates, occur in both the study area such as in entire Central Western Brazil (Funatsu et al., 2012).During the observed period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006), 77% of rainfall precipitations occurred between November through March.The dry period during the austral winter is expressive during June through August, totalizing less than 20 mm of rainfall in average.Maximum monthly precipitations occurred between December and February, whereas maximum runoff showed a delay of one or two months (see Figure 2), having its monthly maxima between February and April (Figure 4a, b).
Visual inspection of the series indicate a decrease of rainfall until the hydrological year 1998/99, followed by a new increase (Figure 4c,  d).The runoff increase in 1999/2000 after the dry year of 1998/99 illustrates the elevated basin memory.One should be aware that local adjustment of the LOESS model to the calendar year (red lines, Figure 4c, d) tend to overestimate runoff at the beginning of the series and to underestimate them at the end.
Both series are characterized by a strong seasonality and have a high serial autocorrelation, which remains significant in a 95% confidence interval, for more than 120 months (Figure 5).Therefore we applied the Kendall test with removal of seasonal influence (MKs) with bootstrapping to minimize the influence of autocorrelation in trend testing (Hipel & Mcleod 2005) (Table 2).
Whereas monthly rainfalls do not have a significant tendency, the MKs test identified a weak positive tendency (tau = 0,176) for monthly runoff, which is not detected.The bias of the bootstrapping statistics, as the difference between the average of 500 resampled sets and the original tau estimate of the MKs test, shows a higher tau of the original series than the average of the resampled series.The absence of a temporal trend in rainfall time series is in accordance with previous climatological studies in the transition zone between the Cerrado and Amazon Forest biomes, which had been little conclusive.
For the Amazon basin, Coelho et al. ( 2013) related a negative tendency after 1976 with a change in the Decadal Pacific Oscillation (Zhang et al., 1997) and Debortoli et al. (2015) identified a reduction in rainfall at about 70% of more than 200 studied gauges, however without detecting significant negative trends in the annual totals.The revision of D'Almeida et al. ( 2007) concluded, that in the Amazon basin no uniform tendencies on major elements of the hydrological cycle can be observed.
At higher latitudes in the mid and southern Cerrado Biome, Collischonn et al. (2001) and Boulanger et al. (2007) identified significant increase in precipitation and runoff.These studies, which evaluated historical data, have been underpinned by model-based scenario projections in the region.Milly et al. (2005), in an ensemble study including simulations with 12 different GMCs, predicted Mato Grosso of being in a transitional zone between Southern regions with precipitation increase and decrease in the Amazon Basin for the period between 2041 and 2060.These projections are in accordance with IPCC (2013), which indicate diminishing precipitations in the Amazon basin and an increase in south-eastern Latin America.
As a non-parametric method, MK tests are not influenced by the magnitude of extreme values.Therefore we analysed as well monthly rainfall and runoff for September, representing the peak of the dry season and of March as the period of maximum runoff and high precipitation.None of these series showed significant serial autocorrelation, neither a significant tendency in MK tests with bootstrapping.
All these results indicate that positive trends in the runoff series are likely to be caused by factors beneath varying rainfall in the basin.

Relationships between runoff and deforestation indicators
Rainfall showed strong positive relationships with yearly runoff in the six subbasins (Figure 6), which is as well positively correlated with the total deforested area (DA).Yearly rates of deforestation (YRD) have an inverse relation with rainfall, which diminished throughout the years, whereas rainfall increased.Thus, if observed independently from the rainfall, the positive relation of YRD with runoff is mimicked.RR increased with DA and therefore in time, governed by the stronger increase of runoff than rainfall.Figure 6.Relationships of yearly rainfall, runoff, runoff ratio and deforestation indicators (DA, YRD) in the six monitored sub-watersheds.Zeilhofer, P., Alcantara, L. H., Fantim-Cruz, I. RR increases when permanency of rainfall in a sub-catchment is reduced, triggered by increasing surface runoff and a reduced sub-surface runoff.
DA shows, with exception of the 17200000 sub-catchment, negative relation with TDA.The two sub-catchments of the Upper Teles Pires (17200000, 17210000) have lower DR values, indicating a low memory, with quicker processes in runoff generation in the steeper declines of the headwaters.
These three independent variables have positive relationships (B coefficients) with the annual runoff.If linear models are adjusted against the years, annual discharge increased in about 13% (~1% p.a.), whereas rainfall decreased in about 2%, hardening the hypothesis that runoff increase was not caused by an alteration of the pluviometric regime.
Our results are hereby in accordance with findings from Coe et al. (2009) in the Tocantins River basin, which was deforested in more than 50% and which showed a similar runoff increase of 1% p.a., this as well without an increase in rainfall.Based on discharge simulations, Coe et al. (2011) concluded that about 2/3 of the runoff increase could be attributed to deforestation in the Tocantins basin.Under tropical semi-humid climates with strong seasonality, an evapotranspiration increase can be expected, if precipitation increases (Costa et al., 2003).Such a pattern was confirmed by Oliveira et al. (2014), who detected an increase of evapotranspiration in the entire Brazilian Cerrado region.Thus, both of these studies underpin our hypothesis that observed runoff increase in the Teles Pires watershed occurred due to deforestation and were not caused by an eventual decrease of evapotranspiration.
Table 3. Likelihood Ratio Chi-square (LRC), B coefficient (B) and Wald χ 2 statistics with the significance levels (*: p<0,05; **: p<0,01) of the independent variables: Rainfall in mm (P) with one month delay, Deforested area in % (DA), Yearly Rate of Deforestation in % (YRD) for GLMs with the annual runoff and runoff ratio as dependants.We denote a scale influence of predictive power of the two evaluated deforestation-related explanatory variables.DA e YRD are significant predictors in the three sub-watersheds smaller than 14.000 km 2 (17200000, 17210000 and 17230000), but loose predictive power in the three larger subwatersheds, which have contribution areas greater than 34,000 km 2 (17280000, 17300000 and 17340000).In the three smaller sub-watersheds the Wald χ 2 statistics of the GLM of annual runoff 3.91 e 18.81 for DA as independent variable, this versus 0.46 e 3.41 in the three larger sub-watersheds (Table 3).For YRD, the Wald χ 2 statistics range between 6.90 and 15.55 in the smaller versus 0.08 and 3.80 in the larger sub-catchments, respectively.A similar decrease of the predictive power of these two groups of sub-watersheds is noted as well in the GLM for RR as dependant variable.Bloeschl et al. (2007) showed that the influence of LUCC decreases with an increase of watershed extension, however with modulations according to other environmental characteristics.Furthermore, they depicted that fine-scaled studies in small paired watersheds, evidenced that runoff increases after deforestation, but that such results cannot be extrapolated to other scales, where the relative importance of runoff generation processes change.Our study in a watershed scale showed, that discharge sensitivity to deforestation decreases exponentially with an increase of contribution area (Figure 7).Both, the Wald χ 2 statistics for DA (Figure 7a) and YRD (Figure 7b) as predictors for annual discharge and RR diminished with increasing watershed size, showing the progressive reduction of the positive relationship between deforestation and runoff.

Conclusions
During the evaluated 15-year period, the upper Teles Pires River watershed has experimented deforestation taxes of about 2.5% p.a. in average, reaching 56.2% of deforested areas in 2007.
During this period runoff showed a significant increase of about 1% p.a., while rainfall increased in less than 0.5%, without showing a significant trend.Annual runoff and the runoff ratio increased with increases in the deforestation indicators DA and YRD.
Magnitude of deforestation impacts on runoff were hereby a function of the watershed size, which are at least six time higher, if the smallest sub-catchment is compared with the largest one.
Therefore we recommend that increased efforts should be undertaken by authorities to monitor and foster sustainable management initiatives in smaller watersheds, where negative impacts of deforestation on runoff are expected to be more expressive than in the great, already monitored hydrological systems.

Figure 1 .
Figure 1.Rainfall stations, river gauges and their contribution areas in the Upper Teles Pires River watershed.Average annual precipitation was interpolated by splines without tension.

Figure 3 .
Figure 3. Deforestation dynamics in the upper Teles Pires watershed between 1992 and 2007, based on the monitoring of FEMA/SEMA.a) Spatio-temporal dynamics and b) deforestation time series in the six gauged contribution areas.

Figure 4 .
Figure 4. Variations of monthly rainfall (a) in the upper Teles Pires watershed and corresponding monthly specific discharge (b) at the outlet (Indeco -17430000).Smoothed curves in the P and Q series (c, d) with a 99% confidence interval were obtained by local weighted regressions (LOESS).

Figure 5 .
Figure 5. Temporal autocorrelation function (ACF) of monthly time series (lag) of rainfall (a) in the upper Teles Pires watershed and of runoff (b) at the outlet (Indeco -17430000).

Figure 7 .
Figure 7. Relationships between the Wald Chi-square of DA (a) and YRD (b) and sub-watershed size in GLMs with annual runoff and RR as dependent variables.

Table 1 .
River discharge gauges with its ANA (National Water Agency) code, discharge, specific discharge and physiographic characteristics.TP: Teles Pires, V: Rio Verde.