On the ability of large-scale hydrological models to simulate land use and land cover change impacts in Amazonian basins

ABSTRACT

The MHD-INPE model was applied in the Ji-Parana Basin, a 30 000 km² catchment located in the southwest of the Amazon Basin which has lost more than 50% of its forest since the 1980s, to simulate land use and land cover change impacts on runoff generation process and how they are related to basin topography. Simulation results agree with observational studies in the sense that fast response processes are significant in sub-basins with steep slopes while in basins with gentle topography, the impacts are most visible in slow-response hydrological processes. On the other hand, the model is not able to capture the dependence of LUCC impacts on spatial scales. These discrepancies are probably associated with limitations in the spatial representation of heterogeneities within the model, which become more relevant at larger scales. We also tested the hypothesis that secondary forest growth should be able to compensate the decrease in evapotranspiration due to forest–cropland or forest–grassland conversion at a regional scale. Results showed that despite the small fraction of secondary forest estimated on the basin, the higher evapotranspiration efficiency of this type of forest counterbalances a large fraction of the LUCC impacts on evapotranspiration. This result suggests that enhanced transpiration due to secondary forest could explain, at least in part, the lack of clear LUCC signals in discharge series at larger scales.

EDITOR D. Koutsoyiannis; ASSOCIATE EDITOR T. Wagener

KEYWORDS:

• land use and land cover change
• large-scale basin hydrology
• hydrological modelling
• Amazon
• secondary forest
• MHD-INPE

1 Introduction

Several processes interacting at different spatial and temporal scales, along with different intensities, affect the hydrological responses of a basin. The extent of hydrological impacts from land use and land cover change (LUCC) and climate variability (CV) is strongly dependent on soil water flow pathways. Therefore, an adequate representation of the runoff generation mechanisms is crucial to assess the hydrological impacts of LUCC and CV on a basin scale.

At a small scale, usually up to 100 km² (Liebscher 1993), impacts on hydrological response are controlled by geo-morphological characteristics of the basin and its channel network. Meanwhile, at a large scale, greater than 10⁴ km² (Liebscher 1993), these impacts are highly controlled by spatial heterogeneity of bio-geophysical characteristics such as climate, vegetation, geology, geomorphology and the landscape fragmentation due to anthropic activity. As a result, changes in hydrological response in the basin from LUCC may be delayed, and its effects do not reach stationary conditions which avoid the determination of stationary precipitation–discharge relationships (Bruijnzeel 2004).

LUCC affects the soil structure of a basin. However, the manifestation of these effects varies depending on the environmental characteristics of the basin. Also, conditions of stability in the hydraulic characteristics of the soil modify LUCC impacts. Previous studies suggest that, in most Amazonian soils, changes in hydrologic characteristics driven by LUCC do not strongly affect the surface flow generation mechanism in a basin (Tomasella and Hodnett 1996, Elsenbeer et al. 1999, Zimmermann et al. 2006). However, LUCC does incrementally impact the volume and the frequency of overland flow events as the area of contribution is expanded (Germer et al. 2010).

Anthropogenic induced landscape fragmentation in large-scale Amazon basins affects evapotranspiration rates on both local and regional scales. Sap flow measurements in forest patches have shown the influence of tree clearing on evapotranspiration rates (Giambelluca et al. 2003), where heat diffusivity is likely to enhance tree transpiration rates near the fragment edge. Also secondary forests, which develop in deforested abandoned areas and are becoming the main type of vegetation in tropical regions (Uhl et al. 1988, Moran et al. 1994), show evapotranspiration rates exceeding those of pristine forests (Shuttleworth 1988, Giambelluca 2002, Von Randow et al. 2010). Hydrological characteristics of secondary forests are highly variable in time (Brown and Lugo 1990, Malmer 1992, Giambelluca et al. 1997). Secondary forest albedo is initially high and then it diminishes reaching pristine forest values after 30 years. In addition, radiative efficiency increases with time (Hölscher et al. 1997, Giambelluca et al. 1997, 1999, Giambelluca 2002). Giambelluca (2002) and Giambelluca et al. (2003) concluded that evaporation rates greater than available radiant energy would be due to an input of advected energy from the surrounding clearings. This implies that fragmented forests and secondary forest growth should be able to compensate for the decrease in evapotranspiration due to forest–cropland or forest–grassland conversion at a regional scale (Giambelluca et al. 2003).

The impacts of LUCC on catchment hydrological response have attracted attention and stirred heated debates among hydrologists for decades, in particular the processes of deforestation and afforestation. A detailed and comprehensive review of this issue has been addressed by Andréassian (2004), Bruijnzeel (2004) and Bonell (2010), among others. Many authors have investigated the effects of LUCC on large-scale basin hydrology finding contradictory results. While several studies conclude that there are no significant changes in the hydrological response of large-scale basins due to LUCC (Buttle and Metcalfe 2000, Wilk et al. 2001, Robinson et al. 2003), other investigations show important effects on discharge (Matheussen et al. 2000, Costa et al. 2003, Siriwardena et al. 2006, Lin and Wei 2008). Results show that the hydrological effects of LUCC in one particular basin cannot be readily extrapolated to another basin, since the effects are scale and time dependent and related to the relief of the drainage area (Rodriguez et al. 2010).

Hydrological models have been used in LUCC studies in order to understand the different mechanisms involved in this process at large-scale basins (Lørup et al. 1998, Vanshaar et al. 2002, Bathurst et al. 2004, Li et al. 2007, Coe et al. 2009, Mueller et al. 2009). In these studies, different land use scenarios, including disturbed and pristine vegetation, are parameterized within hydrological models. Usually, these scenarios are based on pre and post LUCC conditions or on hypothetical future scenarios, such as the complete replacement of pristine vegetation. In general, the results from these simulations represent a state of equilibrium and suggest that significant changes in the hydrological response of a basin are controlled by thresholds directly related to the size of the area affected by LUCC and are associated with changes in evaporative fluxes (Coe et al. 2009).

Different approaches have been proposed to account for LUCC impacts using numerical models, as an alternative to the paired catchment approach, that are applicable to large scales. The use of multi-model ensembles has been proposed as a solution to the issue of dealing with the objective selection of a most appropriate model structure (Breuer et al. 2009, Plesca et al. 2012). Another proposed approach is to calibrate the parameters of a hydrological model using discharge data for pre and post LUCC periods, and then analyse if the changes in the calibrated parameters are significant (Seibert and Mcdonnell 2010, Birkel et al. 2012, Gebrehiwot et al. 2013). Using this approach, Seibert and Mcdonnell (2010) successfully carried out analyses of hydrological change associated to LUCC using a simple conceptual model combined with uncertainty analysis of the calibrated parameters. Limitations of this approach are associated with the availability of representative data pre and post LUCC (Seibert and Mcdonnell 2010). In addition, the model should be able to represent processes associated with LUCC and to integrate those processes across spatial scales (Bonell 2010).

Model responses to LUCC depend on the structure and the parameterizations used in the model. There are two basic methodologies adopted to define the structure of a hydrological model: downward and upward approaches (Klemeš 1983, Dooge 1986). Following Klemeš (1983), the upwards approach is defined as the route that “attempts to combine, by mathematical synthesis, the empirical facts and theoretical knowledge available at a lower level of scale, into theories capable of predicting events to be expected at a higher, in our case hydrological, level”; meanwhile the downwards approach is defined as the route that “starts with trying to find a distinct conceptual node directly at the level of interest (or higher) and then looks for the steps that could have led to it from a lower level”. The upward approach is more appropriate for identifying causal relationships, although the results are highly affected by assumptions used in the development of the model. In addition, model structure and parameter value definitions are strongly affected by scale issues and their inter-relationships. The downward approach is more appropriate for studying the effects of LUCC, but it is more difficult to identify causal relationships (Dooge 1997, Sivapalan 2005).

In this study, we use previous analyses of LUCC impacts in observed discharges in the Ji-Paraná Basin (Rodriguez et al. 2010) to test a downward philosophy model to simulate LUCC impacts on a large-scale basin in the Amazon region. We evaluate the approach’s capacity to accurately represent transient impacts in runoff generation processes at the subgrid-scale and the partition of the LUCC’s impacts between surface and baseflow as described by Rodriguez et al. (2010). We also investigate the ability of the numerical scheme to integrate those processes across time and spatial scales comparing observed and simulated hydrographs in different Ji-Paraná sub-basins. We analyse the model’s skill in simulating observed trends in the flow duration curve along with land use history, specifically forest conversion to pasture scenario including or excluding regeneration (secondary growth). The study area is described in the following section. Model implementation, data used and a description of the experiments are then presented in Section 3. Sections 4 and 5 describe the results of the experiments and offer some discussion. Finally, we present our conclusions.

2 The Ji-Paranâ basin

The Ji-Paraná basin is located in the southwest of the Brazilian Amazon. The drainage area above the Ji-Paraná station (10.87°S–61.94°W) is about 33 000 km² (Fig. 1). The basin is located within the Amazonian Arc of Deforestation and has been losing pristine forest since 1980s (Fig. 2). Agricultural expansion has been the main driver of deforestation in the basin which began in the 1970s (Machado 1998, Alves et al. 2003).

Figure 1. Ji-Paraná Basin geographic location, topography and gauging stations.

Figure 2. Annual deforestation areas as a percentage of total sub-basin area.

A previous study of LUCC in the Ji-Paraná basin concludes that LUCC impact on basin hydrologic response is scale dependent and time delayed (Rodriguez et al. 2010, hereafter RTL2010). This study identified that changes in hydrologic response from LUCC impacts is controlled by physical characteristics of the basin, particularly topography. Changes in the fast hydrologic response, associated with LUCC, were detected in the western branch of the basin, where the topography is the steepest. In the eastern part of the basin, dominated by plateaus, the impacts were mainly found in the slow hydrologic response. Regarding precipitation variability, analyses carried out by RTL2010, Espinoza Villar et al. (2009) and Marengo (2009) did not find any evidence of change in the observed precipitation over this area that could be associated with land use and land cover change. Climate large-scale variability largely drives changes of precipitation at interdecadal scales, and explains the variability observed in hydrological series (Espinoza Villar et al. 2009, Marengo 2009), which has also been shown by Rodriguez et al. (2014) using regional atmospheric simulations.

3 Data and methodology

3.1 Data sets

Daily discharge data from seven stations in the basin were obtained from the Brazilian Water Agency (ANA) (Fig. 1). Gridded hydrometeorological data were obtained from the Center for Weather Forecast and Climate Research (CPTEC) at the Brazilian Institute for Space Research (INPE). Both discharge and hydrometeorological data were thoroughly cleaned to ensure quality datasets (ANEEL 1982, ANA 2006, Rodriguez et al. 2010).

Soil information was extracted from soil surveys of EMBRAPA (1980). Soil profile information, namely texture, was used to estimate hydrological characteristics using pedo-transfer functions (Tomasella et al. 2003). Topographic data for hydrological modelling was obtained from the Shuttle Radar Topography Mission, SRTM (Farr et al. 2007). This was available at 90 m resolution.

Maps of land use classification were obtained from Linhares (2005), for the years 1978, 1980 and 1983–2001, and from the programme of deforestation monitoring of Brasil, PRODES (INPE 2004), available at http://www.obt.inpe.br/prodes/index.php, for the years 2002–2006. This information was previously used by RTL2010 to analyse LUCC impacts in the Ji-Parana basin. PRODES and Linhares (2005) data identify deforested areas in the Brazilian Amazon for each year using LANDSAT images during the dry season to minimize the influence of cloudiness. Annual rates of deforestation are established from August 1st of each year. Also, in PRODES and in Linhares (2005) classification, areas identified as “deforestation” in a particular year, are considered to remain deforested in subsequent years, i.e. secondary regrowth is neglected.

To study the potential effects of secondary growth, fractions of secondary forest (P(SF)) were estimated using the Neef et al. (2006) approach, based on information provided about the fraction of remaining pristine forest (P(PF)). The Neef et al. (2006) statistical model uses a quadratic relationship to estimate fractions of secondary forest (equation (1)) and a negative exponential function to estimate the average age of secondary growth (t_SF) in years (equation (2)) from the remaining pristine forest.

(1)

(2)

Neef et al. (2006) developed this simple statistical model using information from remaining pristine forests. The statistical fit was found using data from previously published works, where secondary vegetation was calculated based on satellite measurements and in situ validations. This data included information obtained over the Ji-Paraná catchment area by Nelson et al. (2000), Roberts et al. (2002), Alves et al. (2003) and Neef et al. (2006).

In our study, we are testing the hypothesis that an enhanced evapotranspiration could reduce the expected impacts of LUCC in discharges. Specifically, we are looking for potential effects that a more water efficient secondary-growing vegetation could have on the hydrology of a large-scale basin that has experienced deforestation. Since there is no secondary vegetation series available for the entire basin, we used the statistical model developed by Neef et al. (2006) to develop hypothetical scenarios for use in hydrological modelling.

Secondary growth is likely to occur on abandoned deforested areas, so for this study, we assumed that the increase of secondary forest areas occur at the expense of the reduction of non-forested areas. Three groups of secondary forest were considered based on regrowth age: less than 5 years old, between 5 and 15 years old and more than 15 years old.

3.2 Hydrological modelling

The grid-based distributed hydrological model of the National Institute for Space Research, MHD-INPE (Rodriguez 2011) was developed from the large-scale hydrological model, MGB-IPH (Collischonn et al. 2007). In the MHD-INPE, the runoff generation processes are represented by a combination of the Xinanjiang Model approach (Zhao 1992, Zhao and Liu 1995) and the TopModel formulation (Beven and Kirkby 1979).

For large-scale basin models, in which each grid point includes high landscape heterogeneity, it is reasonable to assume that the spatial variability of soil depths cannot be explained solely by topography, as in TopModel. The inclusion of additional model parameters is justified both by constrained vertical flows in shallow soils and by TopModel’s limitation in representing water table depth in deep aquifer systems.

The MHD-INPE model includes modules to solve the soil water budget, potential evaporation, transpiration, surface, sub-surface and baseflows in each grid cell. It also includes a channel routing scheme. It is a regular grid distributed model, which is fixed at 5 km in this application. Each grid-cell is sub-divided into hydrological response units (HRU) resulting from the combination of land-uses and soil types. The MHD-INPE model solves the water balance for each unit and calculates an area weighted average flux value for each grid-cell. This tile type approach is commonly used in large-scale hydrological models (e.g. VIC; Liang et al. 1994, MGB-IPH; Collischonn et al. 2007).

In this application, the model was configured to consider two soil layers. The surface layer, usually thinner, generates surface flow and quick responses to atmospheric forcing. The deeper soil layer is used to solve baseflows and it collaborates to keep evapotranspiration rates. Root water uptake from the soil is solved considering Jarvis’ (1989) root distribution approach. Canopy interception is simulated using the Gash et al. (1995) methodology. Simple linear reservoirs are used for flow routing into each cell. Flow exchange between cells is solved exclusively through a channel network using the Muskingum-Cunge algorithm (Cunge 1969).

The MHD-INPE formulation assumes the water table follows the local topographic surface and a non-linear variation of transmissivity in depth (Beven and Kirkby 1979, Iorgulescu and Musy 1997). By assuming that the groundwater flux is at a maximum when water storage in the soil is also at a maximum, the groundwater flux can then be obtained from equation (3):

(3)

where T_sub is the transmissivity when the water table is at the surface; tan β is the slope of the local topography; μ is a parameter that defines the shape of the relationship between transmissivity and soil depth; ξ is the drainable porosity; S_max is the maximum capacity of water storage in the cell and λ_μ is the average contribution area per contour unit, assuming a potential profile for transmissivity (Iorgulescu and Musy 1997). The topographic parameter, λ_μ, is solved independently for each regular grid-cell.

S_max is solved for the bottom soil layer in the model considering its depth, D₃, and soil porosity, ϕ:

(4)

Approaches commonly used in large-scale models assume that any volume of water added in the cell is instantly distributed and the water table is spatially uniform at any moment. On the other hand, the MHD-INPE approach uses the assumption of quasi-permanent state of TopModel to describe the distribution of the water table in the cell.

The percentage of saturated area in the cell, A_sat, is solved through the histogram of the upstream contributing area of the grid–cell. The saturation overland flow, Q_sup, from the saturated area is generated by the precipitation falling on the saturated area:

(5)

MHD-INPE’s structure has an upper soil layer, with uniform depth D₁, which drives the lag-time between infiltration and deep recharge. The sub-surface flux from this layer is solved with equation (6):

(6)

where K_SS is the saturated hydraulic conductivity in the upper soil layer; SS_max is the maximum water storage in this layer; η is the Brooks-Corey parameter; α includes the soil anisotropy in the equation; SS_t is the average water storage in the upper layer and λ_η is the average value of the contribution area per contour unit raised to the 1/η potency.

The drainage from upper to bottom soil layer is solved by equation (7), where λ represents the upward contribution per unit contour:

(7)

Table 1 shows the model parameters considered for calibration. Prescript values for vegetation and soil parameters were obtained from literature (Jones 1992, Shuttleworth 1993, Culf et al. 1996, Roberts et al. 1996, Mc William et al. 1996, Ubarana 1996, Wright et al. 1996, Giambelluca 2002, Aragao et al. 2005).

Table 1. List of calibrated model parameters.

Parameter	Unit	Physical meaning
D1	m	Thickness of the upper layer
D3	m	Thickness of the bottom layer
Kss	m d^-1	Saturated hydraulic conductivity
Tsub	m^2 d^-1	Maximum transmissivity of the bottom layer
µ	–	Parameter that represents the decay of transmissivity with the thickness of the saturated zone
Csup	d^-1	Routing water storage parameter for surface flows
Csub	d^-1	Routing water storage parameter for subsurface flows

We use the Paz and Collischonn (2007, 2008) approach to correct the SRTM data and to obtain the drainage network. The methodology developed by Paz and Collischonn (2007, 2008) uses globally available DEM data with higher resolution to derive river geomorphologic characteristics at the model grid-cell scale by upscaling procedures, preserving the total stream length. The Quinn et al. (1991) methodology was used to calculate the topographic index distribution in each of the model’s cells. Quinn et al. (1991) develops a multiple flow direction algorithm to determine the distribution of cumulative area and topographic index, which distributes the accumulated upslope area for each cell among all downslope directions. In order to obtain similar HRUs, we group soil types into four groups within the basin following similarities in hydrological characteristics (Doyle et al. 2013). We then combine these with land use classifications in a 100 m × 100 m resolution, resulting in eight HRUs for the basin.

Following RTL2010, we carried out statistical analyses of the simulated daily discharges for each hydrologic year, starting in October and finishing in September. For each flow duration curve (FDC), we obtain the average discharge (AVG); median (MDN); lower (LWQ) and upper (UPQ) quartiles; second (P02), 10th (P10) and 90th (P90) percentiles; variance (VAR), skewness (SKW) and kurtosis (KUR). Also, the signatures of the FDC described by Yilmaz et al. (2008) were obtained: the mid-segment slope of the FDC (SLP, equation (8)), the volume of water belonging to the high-flow segment of the FDC (VMX, equation (9)) and the total volume belonging to the low-flow segment of the FDC (SLP, equation (10)). The SLP signature was defined considering the 20% exceedance probability discharge (Q_m1) and 70% exceedance probability discharge (Q_m2), where m1 and m2 are the respective points in the FDC. VMX was evaluated for discharges with exceedance probability lower than 2% (m3), and VMN for discharges with exceedance probability larger than 70% (m2). In equations (8)–(10), Δf is the frequency increment in the FDC and, in equation (10), Q_i is the ith exceedance probability discharge.

(8)

(9)

(10)

3.3 Numerical experiments

Because historical discharge data in the basin began concomitantly with agriculture expansion, it is not possible to separate the simulation period into pre and post LUCC. Because of this, the experiments were designed to include land use and land changes history in the simulation period. Three experiments were developed with the MHD-INPE model: EPAST, EFOR and ECAPO. In the EPAST experiment, LUCC between 1982 and 2006 are obtained from annual vegetation surveys. Because of the lack of land use and land cover information and the small variation in landscape conditions between 1980 and 1983 in the basin (Rodriguez et al. 2010), we assume the condition for 1982 to be as in 1980. This experiment was used to calibrate the hydrological model. The EFOR experiment uses vegetation coverage fixed to 1978 conditions, when non-forested areas in the basin were less than 5%. Finally, in the ECAPO experiment, secondary forest percentages obtained from the Neef et al. (2006) model are included.

Secondary forest evapotranspiration in the ECAPO experiment was parameterized in the MHD-INPE model to reproduce the behaviour described by Giambelluca (2002), with a higher evaporative fraction than primary forest. In young secondary forest, albedo is higher than in pristine forest. Albedo decreases with the age of secondary forest, until it reaches values similar to pristine forest after 30 years. In order to reproduce this behaviour in the model, albedo values for secondary forest were set as explained in Giambelluca (2002), while a minimum surface resistance was used to fit the evaporative fraction behaviour.

Calibration was performed through a manual approach, using both performance coefficients and visual inspection of results. To assess the goodness of fit between observations and simulations, five performance coefficients were used in this study: the Nash-Sutcliffe model efficiency, NSE (equation (11)) (Nash and Sutcliffe 1970); the Kling-Gupta coefficient, KGE (equation (12)) (introduced by Gupta et al. 2009); the determination coefficient using both discharges, R² (equation (13)) and the logarithm of discharges, LR² (equation (14)); and the bias, expressed as differences between simulated and observed volumes, BIAS (equation (15)). The NSE and LR² coefficients were considered for the calibration process, while the other coefficients were used to carry out a more complete analysis of performance. We applied a simple approximation of the Pareto solution considering both coefficients (Vrugt et al. 2003). Since the goal of this work is to analyse model capacity to reproduce the variability observed in the main features of the discharge duration curve associated with LUCC (Rodriguez et al. 2010), we evaluated the model performance for the entire simulated period.

(11)

(12)

(13)

(14)

(15)

where QS_i is the simulated discharge at time ; QO_i is the observed discharge at time ; r is the linear correlation coefficient; σ_S is the standard deviation of simulated discharges; σ₀ is the standard deviation of observed discharges; μ_S is the mean of simulated discharges and μ₀ is the mean of observed discharges.

The NSE, as defined here, is a normalized measure, varying between −∞ and 1, which compares model performance with the performance of a simple base-line model that uses a constant mean observed value as a discharge estimator. A negative value of NSE means the model performance is worse than the performance obtained using the base-line model, and values greater than 0 are associated with an acceptable model performance (Moriasi et al. 2007). On the other hand, maximizing NSE implies an underestimation of variability (Gupta et al. 2009).

Gupta et al. (2009) propose a performance criterion, termed the KGE, which improves the bias and variability measure. The KGE compares the Euclidean distance from the ideal point of three components related to correlation. The KGE takes into consideration bias, expressed as the ratio between observed and simulated means, and variability, expressed as the ratio between observed and simulated standard deviations (Gupta et al. 2009).

R² is a measure of how much of the observed dispersion is explained by the model. R² varies between 0 and 1, and, in general, values greater than 0.5 are considered acceptable (Moriasi et al. 2007). High values of R² mean that the model is able to assess the observed dispersion (Krause et al. 2005). To take into account low discharges, the coefficient of determination is calculated for the logarithm of discharges, LR². This provides a balance between criteria that penalise errors at high flows and others that are less sensitive to a small number of large errors at high flows. BIAS is a measure of the fit between observed and simulated volumes. It has been included here as a reference value, once LUCC’s impact will be also described over the annual total discharge.

4 Results

4.1 Model performance

Model performance should be evaluated within the context of the main goals for which its use is destined (Schaefli and Gupta 2007). In this work we are interested in assessing the model’s ability to simulate observed LUCC’s impact rather than in its capacity to predict discharge for the basin.

First, we looked at the representation of variability and the explained dispersion in the model for the entire simulation period. Performance coefficients (Table 2) indicate that the model’s performance in simulating historical discharges in the basins is satisfactory (Moriasi et al. 2007). The first two years of the simulation period (1982–1983) are neglected since they are treated as part of the model warm-up. The lowest NSE value is found at SB3, however it is still greater than 0.5. R² and LR² show that the model is able to explain more than 65% of the observed dispersion. KGE also indicates a good performance in capturing variability. The error in simulated volume, represented by BIAS, is lower than 15%, an error usually acceptable for the stage-discharge curve. The MHD-INPE model produces an adequate simulation of the recession curve and baseflow contribution (Fig. 3); though it has difficulty capturing the early recovery in discharges at the end of the dry season and in matching storm event peaks (Fig. 3). These limitations might be related to TopModel’s assumptions such as constant recharge and the quasi-steady state assumption. In areas where the water table is close to the surface, it is reasonable to expect a quick recharge during intense storms, while recharge would be expected to be delayed in areas where the water table is deeper. Areas with a shallow water table are commonly adjacent to streams. Therefore, an intense recharge in these areas could result in a fast response in channels, contributing to high peaks and to the early recovery of the hydrograph.

Table 2. Performance coefficients for daily discharges (NSE, KGE, R² and LR²) and volumes (BIAS) for six sub-basins of the Ji-Parana basin.

	NSE	KGE	R^2	LR^2	BIAS (%)
SB1	0.63	0.80	0.71	0.73	−8.98
SB2	0.73	0.85	0.87	0.78	−10.75
SB3	0.54	0.76	0.66	0.72	1.24
SB4	0.67	0.78	0.78	0.78	−7.75
SB5	0.77	0.88	0.85	0.85	−3.48
SB6	0.80	0.89	0.86	0.86	−2.74
SB7	0.83	0.89	0.80	0.90	3.43

Figure 3. Hydrographs for observed (black line) and simulated (red line) discharges at (a) SB2, (b) SB5 and (c) SB7.

There is a gap in the quality of the meteorological data before and after 1996 in the basin (Rodriguez et al. 2010). This gap affects the model performance when evaluating different time periods. During 1984–1995, the model performance is worse than during the period 1996–2005. The lowest NSE during 1984–1995 is 0.52, in the SB3 sub-basin; while during 1996–2005, it is 0.64, in the SB1 sub-basin. It should be noted that, in the SB3 sub-basin, discharge records are poor during 1984–1998. Such behaviour could affect the model’s skill in simulating trends observed in the series of hydrological parameters as described by RLT2010.

4.2 Impacts of forest to grassland changes

The EPAST experiment shows an increment in annual discharge in the basin when deforestation is considered in the simulations (Figs 2 and 4). Increments in discharge are associated with decreases in evapotranspiration driven by forest-to-grassland changes (Table 3). Since pasture evapotranspiration is lower than forest evapotranspiration, discharges are higher under EPAST conditions, and discharge increases proportionally as deforested area also increases. The Mann-Kendall test (Sneyers 1975) identifies statistically significant positive trends in annual discharge differences between EPAST and EFOR at a 95% confidence level. The Kendall non-parametric test (Kendall 1975) shows a statistically significant correlation between these differences and the percentage of deforestation in the basin (not shown).

Table 3. Changes in evapotranspiration (ET) and annual total discharge (VOL) between experiments considering pristine conditions (EFOR) and forest-to-grass change (EPAST).

	SB1	SB2	SB3	SB4	SB5	SB6	SB7
ETEPAST—ETEFOR (mm)	21.13	28.71	39.82	73.85	75.10	57.93	102.31
VOLEPAST—VOLEFOR (mm)	−21.61	−28.93	−40.71	−74.32	−76.52	−58.64	−103.22

Figure 4. Annual total discharge differences between experiments considering pristine conditions (EFOR) and forest-to-grass change (EPAST). Values shown as a percentage of annual total discharge under pristine conditions ([EPAST-EFOR]/EFOR*100).

By simulating LUCC in the basin, the model reveals that surface and baseflows are affected differently in the eastern and western branches of the upper basin. While in sub-basins SB1 and SB2, in the eastern upper basin, the impacts are greater in baseflow than in surface flow, in sub-basins SB4 and SB5, on the western side, the impacts are higher on surface flux (Fig. 5).

Figure 5. Contribution of baseflow and surface flow impacts for annual total discharge differences under forest-to-grass change in (a) eastern SB2 and (b) western SB5 sub-basins. Values shown as a percentage of total difference explained by difference in each flow.

On the other hand, simulated discharges do not show the same trends in statistical parameters as shown by RTL2010 for the observed discharges during the historical period (Table 4), In the RTL2010 study negative trends in statistical parameters related to the lowest discharges and positive trends related to the largest discharges in the smaller sub-basin. Looking at the same periods as those analysed by RTL2010, our simulated series show statistically significant positive trends mainly in parameters associated with average to low discharges (Table 4), which does not agree with RTL2010’s results.

Table 4. Statistically significant trends in model parameters from observed and simulated daily discharges series. “+” denotes positive trends and “−” denotes negative trends.

	Observations							EPAST (Simulations)
	SB1	SB2	SB3	SB4	SB5	SB6	SB7	SB1	SB2	SB3	SB4	SB5	SB6	SB7
	91–04	90–99	91–95	91–04	90–04	95–04	93–01	91–04	90–99	91–95	91–04	90–04	95–04	93–01
AVG				+							+
MDN	−										+	+
UPQ								+
LWQ	−	−			−						+			+
P90	−	−						+			+			+
P10				+
VAR				+
SKW	+		−							−
KUR	+									−
P02
SLP
VMX	+
VMN	−

4.3 Impacts of secondary forest

Secondary forest areas, simulated with the Neef model (Neef et al. 2006), represent a small percentage of the total basin area (Fig. 6(a)), which is supported by observations in the region (Nelson et al. 2000, Roberts et al. 2002, Alves et al. 2003). Secondary vegetation parameters are set to reproduce the increment in the evaporative fraction described by Giambelluca (2002) (Fig. 6(b) and Table 5). This approach results in evapotranspiration rates higher than those observed in pristine forest (Table 5). The results are in agreement with descriptions for secondary forest transpiration at an experimental site in central Amazon described in Von Randow et al. (2010).

Table 5. Daily averaged evapotranspiration (ET) and evaporative fraction (EF) for different vegetation types.

TIPO	ET (mm/d)	EF
Non-forest	2.70	0.91
Forest	3.39	0.99
Secondary vegetation
<5 years old	3.56	1.06
5–15 years old	3.74	1.10
>15 years old	3.86	1.12

Figure 6. (a). Total secondary vegetation fraction in each sub-basin and (b) evaporative fraction for pasture (grey line), forest (black line), <5-year-old secondary vegetation (c1), 5–15-year-old secondary vegetation (c2) and >15-year-old secondary vegetation (c3).

The model experiment that accounts for secondary growth forest results in a higher annual total discharge than the model experiment that only considers permanent pristine conditions (Fig. 7(a)). On the other hand, when comparing the model experiments for secondary growth forests with forest-to-grassland conversion; secondary growth forest results in lower annual total discharge values (Figs 4 and 7(a)). Also, the influence of secondary vegetation on annual total discharge shows a gradual reduction over time (Fig. 7(b)). This is directly related to the fact that during the simulation period deforested areas expanded much more rapidly than the regeneration of secondary forest growth. The experiment reveals that, in spite of the small areal percentages, secondary forest growth has an important impact on annual total discharge. At the end of the simulation period, when deforestation rates are large, secondary forest growth still reduces the forest-to-grassland impact on annual total discharge by more than 20% (Fig. 7(b)). In the upper basin, at the beginning of the simulated period when deforested areas are small, annual total discharge under secondary forest growth is smaller than under pristine conditions (Fig. 7(b)).

Figure 7. (a). Annual total discharge differences between experiments considering pristine conditions (EFOR) and forest-to-grass change and secondary vegetation growth (ECAPO). Values shown as percentage of annual total discharge under pristine conditions ([ECAPO-EFOR]/EFOR*100) for each sub-basin. (b) Reduction in the forest-to-grass change impacts due to secondary vegetation growth for each sub-basin.

LUCC also affects the extension of variable source contribution areas in the basin (Fig. 8). Under the forest-to-grassland change scenario, in the EPAST experiment, the extension of saturated area remains larger than in the EFOR experiment. When secondary vegetation growth is considered, the impact on the saturated area is smaller than in the EPAST experiment. The variability of the saturated area is smaller in the western sub-basins than in the eastern ones (Fig. 8). These differences affect both runoff generation and soil evaporation, which have opposite effects in annual total discharge.

Figure 8. Exceedance probability of daily saturated area (%) in (a) eastern SB2 and (b) western SB5 sub-basins for experiments under pristine conditions (EFOR), forest-to-grass change (EPAST) and secondary vegetation growth (ECAPO).

5 Discussion

Contrasting with the numerical experiments generally found in literature, which are based in static land-use scenarios, we used the historical series of land use as input in a hydrological model and compared simulated discharges with observations. The objective of this experiment was to determine if the approaches used in large-scale models were able to simulate the scale dependence of LUCC impacts. Time-lagged effect of forest conversion could be an explanation for the lack of signals in hydrological series. In other words, since previous modelling experiments were conducted for equilibrium conditions (Lørup et al. 1998, Vanshaar et al. 2002, Bathurst et al. 2004, Li et al. 2007, Coe et al. 2009, Mueller et al. 2009), we tested the model ability to reproduce the gradual changes of the hydrological response associated to land conversion.

To verify the performance of the hydrological model, we applied a variety of performance coefficients which capture different behaviours of the time series variability. Performance coefficients used in this study are standard values applied in the use of hydrological models, included in experiments addressing LUCC impacts. All of them suggest that the performance of the model is acceptable according to the standard thresholds reported in literature (Moriasi et al. 2007). The fitting of the MHD-INPE simulations to historical discharges for each sub-basin is acceptable for the goal of this study and similar to the fit reached by other studies addressing numerical simulations of LUCC impacts (Breuer et al. 2009, Li et al. 2009, Mao and Cherkauer 2009, Van Roosmalen et al. 2009, Elfert and Bormann 2010, Ghaffari et al. 2010).

On the other hand, MHD-INPE results do not show the trends reported by the RTL2010 study in series of hydrological parameters obtained from observed duration curves. However, meteorological data quality could have affected the ability of the model to simulate the observed trends in the characteristics of the discharge duration curves.

MHD-INPE simulation of forest-to-grassland changes in the Ji-Paraná basin results in increases in annual total discharge due to decreases in evapotranspiration. Annual total discharge differences have a general positive trend, increasing with deforestation although they show interannual variability likely to be associated with climate conditions over the basin. Such a trend is observed in the results from all sub-basins across all scales. The model’s results show nearly linear impact increases with basin scale, while analysis based on observed data indicates that LUCC impacts are diluted at larger scales (Rodriguez et al. 2010). These limitations for reproducing the observed behaviour in the Ji-Paraná catchment are also verified in the VIC model (Santiago 2005) and in the MGB—IPH model (Bayer and Collischonn 2013).

These discrepancies are probably associated with limitations in the spatial representation of heterogeneities within hydrological models, which become more relevant at larger scales, enhancing transient processes associated to the large variety of land-use transition, changes in soil properties, vegetation growing stage, influence of groundwater, among others. For instance, land use sub-grid variability is identified as a percentage of each land use category within the grid-cell. The approach does not take into consideration the location of each land use within the grid-cell, nor the interaction between land uses. Effective fluxes for each grid-cell are calculated as a weighted average of the total fluxes from each of the HRUs that correspond to a particular grid-cell. The juxtaposition of patches with different hydrological behaviour could have a buffering effect on the impacts of LUCC. For example, runoff may be generated in a highly impervious deforested patch but then could be re-infiltrated when passing through a forested patch with higher infiltration capacity (Ziegler et al. 2007). Giambelluca et al. (2003) found that this juxtaposition could enhance transpiration rates at the edge of forest patches. Giambelluca (2002) also found that the growth of secondary vegetation in previously deforested and then abandoned areas could increase transpiration in the basin. These processes are not included in the model approach; nor are they considered in most large-scale hydrological models. Also, the approaches applied in these models enhance linear scale transference of the LUCC’s impacts. That implies that results of simulated LUCC impacts carried out with these approaches could be overestimated.

We tested the hypothesis that secondary vegetation growth counterbalances expected LUCC impacts. Secondary forest growth diminishes the impact of forest-to-grassland change by triggering higher evapotranspiration rates than primary forest. This could help to explain the lack of signals in the discharge series of the whole basin. Differences in simulated annual total discharge are significantly higher when more than 20% of the basin area is considered to be deforested and when only the forest-to-grassland change is taken into account. This threshold increases to 30% when the effects of secondary vegetation growth are considered; despite the fact that the fraction of the area covered with this vegetation type in the basin is relatively small (Fig. 6(a)). This implies that higher transpiration rates attributed to secondary forest regrowth in abandoned areas counters the increase in annual total discharge due to the loss of forest in the model. This suggests that enhanced transpiration in secondary forest regrowth could be a possible explanation for the lack of impact in discharges from large-scale basins that have partially lost their primary forest. It could also be possible that edge effects on tree transpiration could counterbalance the reduced evapotranspiration rates resulting from forest lost (Giambelluca et al. 2003). RTL2010 shows that edge density increases in the basin during most of the simulation time. The counter effect on LUCC impacts to basin discharge would look similar to those of an increased fraction of secondary vegetation.

In the MHD-INPE’s simulations, LUCC affects differently the surface and baseflow in the eastern and western part of the upper basin. At a cell-model level, these differences in the simulated impacts are explained by differences in the extension of saturated areas. In the eastern sub-basin, where topography is gentler and deforestation rates lower, forest-to-grassland changes result in small differences in saturated area extension, and affect mainly the baseflow; conversely in the western sub-basin, which has a steep topography, forest-to-grassland changes mainly impact fast response generation. These results agree with RTL2010, where observations show that LUCC’s impact is likely to affect lower discharges in the eastern sub-basins and higher discharges in the western sub-basin.

These results show that MHD-INPE is able to simulate the impacts of land use change in runoff generation processes and the dependence of these processes on basin characteristics at the local scale. However, linear integration at the basin scale and misrepresentation of interactions between land use patches at the subgrade may cause the model to lose the ability to capture the scale dependence of LUCC’s impacts.

The experiments carried out do not consider any variation in hydraulic properties of soils after clearance. The potential effects of deforestation on soil properties have been studied in the region by Zimmermann et al. (2006) and Germer et al. (2010). Germer et al. (2010), for instance, concluded “the transformation of undisturbed forest into pasture increased the frequency of occurrence of perched water tables and overland flow by (subsurface) soil compaction and an increase in hydraulic conductivity”.

Germer et al. (2010) studied the effects of deforestation of a large cattle ranch. This land use history differs from the majority of the land uses of the lower Ji-Paraná Basin which is dominated by small-scale properties for subsistence agriculture (Alves et al. 2003). Small properties are cleared using the traditional method of slash and burning (Alves et al. 2003) which causes minimum disturbance in soil properties (Tomasella and Hodnett 1996), while ranches are affected by soil compaction by cattle. Since there is no field information available for all the soil types and varieties of LUCC in the whole basin, the parameterization of soil physics processes cannot be extrapolated at larger scales based on local scale experiments (Gupta et al. 2014).

In addition, considering the high rates of deforestation in the basin and that deforestation rates were more severe in the lower Ji-Paraná basin (Fig. 2), it follows that the impact of land conversion should have been clearly detected at the gauging stations close to the basin mouth. In spite of this, RTL2010 showed that no signal due to LUCC could be detected in the discharge series of the whole Ji-Paraná catchment, although signals are evident in headwater catchments with small drainage area. Therefore, observation suggests that the impacts of subsurface soil compaction observed by Germer et al. (2010) are diluted at larger scales (Rodriguez et al. 2010). The lack of signal at larger scales is consistent with other previous studies for the Amazon Basin: neither Espinoza Villar et al. (2009) nor Marengo (2009) found signals of hydrologic changes that could be associated to deforestation in the analyses of discharges series at several gauging stations in large Amazon sub-basins.

Reductionist paradigm application in large-scale hydrological modelling is based on the use of local scale parameters upscaled to represent processes at basin scale, which results in higher uncertainties of model simulations (e.g. Sivapalan 2005, Gupta et al. 2014). It is important to note that the set of experiments focused on the ability of the model approach to simulate LUCC rather than the uncertainties in calibrated parameters. The input evaporation parameters used in evaporation estimations in the model were independently derived and validated from field evaporation measurements in micrometeorological towers, and they are considered the most realistic representation of this process. Because processes such as the interactions between land uses within the grid-cell, which could explain the lack of the dilution of LUCC impacts with the scale shown by the observations, are not represented by the model formulation, it follows that including parameters uncertainties in the analysis would not provide further answers to model limitations on representing LUCC impacts.

Including processes that enhanced evapotranspiration, such as the growth of secondary vegetation, could help to explain the dilution of the LUCC’s impacts in a large-scale basin with high landscape heterogeneity as suggested by Giambelluca (2002). However, other processes that could also have an important impact in LUCC effects, such as the borders effects in evapotranspiration on forest patches (Giambelluca et al. 2003) or the buffer effect on the overland flow (Zimmermann et al. 2006), are not incorporated in large-scale models.

Since the use of a hydrological model to predict LUCC impacts relies on the correct representation of the relevant processes (Seibert and Mcdonnell 2010), this study suggests that most large-scale hydrological models are unable to represent those impacts in Amazonian basins because of the simplifications adopted in their formulation, which limits the ability to represent the desired processes. Therefore, the main challenge is to incorporate these effects in large-scale models conserving parsimony.

6 Conclusion

We applied the MHD-INPE model, which combines the TopModel and Xinanjiang model philosophies in order to reproduce the impacts of land use and land cover change on runoff generation processes, and their dependence on scale and basin characteristics. Results from model simulations agree with results found by the RTL2010 study that LUCC impacts fast response processes in sub-basins with steep slopes due to an extension of contributing areas, but impacts slow response processes in sub-basins with gentle topography as a result of changes in transpiration losses from deep soil layers.

However, the correct representation of the runoff generation at the local scale (cell-grid model) does not imply the correct representation of processes at larger scales. Moreover, models are able to simulate most of the effects of LUCC at a small scale but fail to upscale such impact at larger scales due to limitation on their formulations. The MHD-INPE model has not been able to capture the scale dependency of the impacts nor the trends in hydrological parameters as shown in the RTL2010 study. The simulation of increasing impacts in conjunction with increasing basin scale is likely to be associated with the misrepresentation of interaction processes between land use patches at the subgrid scale and the linear integration of non-linear hydrological responses through the scale. For instance, how local effects of LUCC are propagated though the drainage network downstream is not fully understood (Pattison and Lane 2012). Most large-scale hydrological models suffer from such limitations.

In this study, we explored if the presence of secondary vegetation, where evaporation rates are higher than primary forest (based on field experiments), could explain the reason for the absence of signals at larger scales, due to the fact that the reduction of evaporation in pasture areas are compensated by secondary forest. Based on modelling results, secondary forest growth moderates the impact of forest-to-grassland conversion on the hydrologic responses of the basin via increasing evapotranspiration losses. This conclusion, based on the assumption that secondary forest growth has a higher evaporative efficiency than pristine vegetation, suggests that processes enhancing transpiration provide a reasonable explanation for the lack of a clear LUCC signal in discharge series at the Ji-Parana station.

For a comprehensive simulation of the transient responses in the hydrology of basins due to LUCC, a better understanding of both the interaction processes between forest and non-forest patches at the local scale and the integrated basin responses at various scales is still necessary. The current approaches in large-scale hydrological models usually predict an equilibrium type response to LUCC and they are unable to represent changing behaviour within a dynamic system. Generating an adequate representation of such behaviour would be useful for identifying critical thresholds, or tipping points, beyond which the system could shift to an alternative configuration, and would provide a scientific base for the detection of early warning signs.

Unfortunately, the development of physically realistic and transposable models, the “holy grail” of hydrology (Gupta et al. 2014), is still the aim of hydrologists. As mentioned by Blöschl et al. (2007) hydrology is a context-dependent discipline. This prevents the “direct” application of the proposed approach to other basins. Therefore, further tests are necessary in other catchments undergoing severe land use changes in the recent past to be able to make a more general assessment of impacts associated to LUCC.

McIntyre et al. (2013) pointed out that “the ability of hydrologists to quantify the impact of rural land use change on the water cycle is however limited and we are not able to provide consistently reliable evidence to support planning and policy decisions. This shortcoming stems mainly from lack of data, but also from lack of modelling methods and tools”. Our results contribute to understanding the limitations of model parameterizations used in large-scale atmospheric and hydrological models. It is clear from our study that current parameterizations used to simulate the LUCC impacts are missing several buffer effects and might be over-predicting the expected changes related to LUCC.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work has been supported by Fundação de Amparo à Pesquisa do Estado de São Paulo – FAPESP process 2008/50436-1.