Modeling freshwater quality scenarios with ecosystem-based adaptation in the headwaters of the Cantareira system, Brazil

Although hydrologic models provide hypothesis testing of complex dynamics occurring at catchments, freshwater quality modeling is still incipient at many subtropical headwaters. In Brazil, a few modeling studies assess freshwater nutrients, limiting policies on hydrologic ecosystem services. This paper aims to compare freshwater quality scenarios under different land-use and land-cover (LULC) change, one of them related to ecosystem-based adaptation (EbA), in Brazilian headwaters. Using the spatially semi-distributed Soil and Water Assessment Tool (SWAT) model, nitrate, total phosphorous (TP) and sediment were modeled in catchments ranging from 7.2 to 1037 km². These headwaters were eligible areas of the Brazilian payment for ecosystem services (PES) projects in the Cantareira water supply system, which had supplied water to 9 million people in the São Paulo metropolitan region (SPMR). We considered SWAT modeling of three LULC scenarios: (i) recent past scenario (S1), with historical LULC in 1990; (ii) current land-use scenario (S2), with LULC for the period 2010–2015 with field validation; and (iii) future land-use scenario with PES (S2 + EbA). This latter scenario proposed forest cover restoration through EbA following the river basin plan by 2035. These three LULC scenarios were tested with a selected record of rainfall and evapotranspiration observed in 2006–2014, with the occurrence of extreme droughts. To assess hydrologic services, we proposed the hydrologic service index (HSI), as a new composite metric comparing water pollution levels (WPL) for reference catchments, related to the grey water footprint (greyWF) and water yield. On the one hand, water quality simulations allowed for the regionalization of greyWF at spatial scales under LULC scenarios. According to the critical threshold, HSI identified areas as less or more sustainable catchments. On the other hand, conservation practices simulated through the S2 + EbA scenario envisaged not only additional and viable best management practices (BMP), but also preventive decision-making at the headwaters of water supply systems.

1 Introduction

Basin plans comprise the main management tool and they plan sustainable use of water resources in both spatial and temporal scales. For sustainable water allocation, river plans are based on accurate data on actual water availability per basin, taking into account water needs for humans, environmental water requirements and the basin's ability to assimilate pollution (Mekonnen et al., 2015). However, adaptive management options such as ecosystem-based adaptation (EbA; see CBD, 2010; BFN/GIZ, 2013) and the water footprint (WF) (Hoekstra and Chapagain, 2008; Hoekstra and Mekonnen, 2012) have rarely been incorporated into Brazilian basin plans. Moreover, integrated qualitative and quantitative simulations and indicators of human appropriation of freshwater resources are seldom used in river plans. The concept of ecosystem-based adaptation is addressed as “using biodiversity and ecosystem services to help people adapt to the adverse effects of climate change”, which was defined by the Convention on Biological Diversity – 10th Conference of the Parties (CoP) (CBD, 2010). Detailed definitions of EbA applied to the Cantareira system's headwaters can be found in Taffarello et al. (2017). The WF is still an environmental indicator used in watershed plans. For example, Spain uses WF as an indicator in basin plans (Hoekstra, 2017; Velázquez et al., 2011; Aldaya et al., 2010). The clean water plan of Vancouver (June 2011) established the reduction of the WF as a sustainable action in its water resources management (MetroVancouver, 2011; Zubrycki et al., 2011). The Colombian government was the first to publish a complete and multi-sectorial evaluation of WF in its territory. Although this study, titled Estudio Nacional del Agua (Colombia, Instituto de Hidrología, Meteorología y Estudios Ambientales, 2014), was not included in the national water management plan, the strategic plan of the Magdalena–Cauca basin incorporates the grey water footprint (greyWF) to assess agriculture pollution (Colombia, 2014). In Brazil, a glossary of terms released by the Brazilian National Water Agency (ANA, 2015) includes the concept of WF to support water resources management.

The WF (Mekonnen and Hoekstra, 2015; Hoekstra et al., 2011) measures both the direct and indirect water use within a river basin. The term water use refers to water withdrawal, as the consumptive use of rainwater (the green water footprint) and of surface and groundwater (the blue water footprint), and water pollution, i.e., the flow of water used to assimilate the pollutant loads (the grey water footprint; see Chapagain et al., 2006). Given that water pollution can be considered a non-consumptive water use, the greyWF is advantageous by quantifying the effects of pollution by flow instead of by concentration, making water demand and availability comparable.

Water footprint assessment comprises four phases: (1) setting goals, (2) accounting, (3) sustainability assessment and (4) response formulation. At the WF response formulation phase, the EbA options, represented by best management practices (BMP) at the catchment scale, could represent a tradeoff on greyWF (Zaffani et al., 2011). That is, BMP adopted in the catchment scale could contribute indirectly to decreasing the level of water pollution. Thus, the EbA would compensate for the greyWF of a certain river basin (Taffarello and Mendiondo, 2011). In the context of water security associated with land-use and land-cover (LULC) change, many existing conflicts over water use could be prevented (Winemiller et al., 2016; Aldaya et al., 2010; Oki and Kanae, 2006). For example, LULC influences water quality, which affects the supporting¹ and regulating² ecosystem services (Mulder et al., 2015; MEA, 2005) and needs to be monitored for adaptive and equitable management on the river basin scale (Taffarello et al., 2016a, b). In spite of discussions regarding the lack of representativeness of data used in early studies with greyWF (Wichelns, 2015; Zhang et al., 2010; Aldaya et al., 2010; Aldaya and Llamas, 2008), we argue that the greyWF method may account for hydrologic services and provide a multidisciplinary, qualitative-quantitative integrated and transparent framework for better water policy decisions. Understanding these catchment-scale ecohydrologic processes requires not only low-frequency sampling, but also automated, in situ high-frequency monitoring (Bieroza et al., 2014; Halliday et al., 2012), as well as using ecohydrologic modeling to protect water quality and quantity. However, freshwater quality modeling associated with EbA, greyWF and LULC is still incipient in many river catchments. In Brazil, approximately only 5 % of modeling studies evaluate nutrients in freshwater (Bressiani et al., 2015), which limits the policies on regulating ecosystem services.

In this research, we propose that the regulating ecosystem services be addressed by the greyWF because it considers the water volume for self-purification of receiving water bodies affected by pollutants (Zhang et al., 2010). The working hypothesis of the paper is related to how conservation practices addressed by EbA impact hydrology and the ecosystem services, such as maintaining, restoring or improving both the water yield and the freshwater quality, using ecohydrologic modeling in different catchment scales. On the other hand, we hypothesized that incentives of EbA policies can affect water yield and water quality through nonlinear tradeoffs, with high spatiotemporal complexity, which can be assessed by modeling, but previously supported by in situ monitoring variables for setup boundary conditions of simulation runs. In these scales, the greyWF can evaluate the changes in the regulating hydrologic services. Among the three water footprint components, in this study we assessed greyWF for nitrate, total phosphorous (TP) and sediments in 20 sub-basins in the headwaters of the Cantareira water supply system. Thus, the aim of this study is to compare freshwater quality scenarios, one of them related to EbA options through BMP, and to assess greyWF under different LULC changes: (S1) historic LULC of 1990, (S2) current LULC for the period 2010–2015, and (S2 + EbA) future LULC based on EbA with S2 as a baseline. This method is addressed using nested catchment experiments (NCE, see Taffarello et al., 2016a, b) at a range of scales from small catchments of 7.7 km² to medium-size basins of 1200 km² at subtropical headwaters responsible for the water supply of the São Paulo metropolitan region (SPMR).

Therefore, this paper consists of four sections. The first section provides a brief description of the context, gap, hypothesis and our research goals. The second section describes the simulation methods used in the watershed scale and the development of three LULC scenarios. We then propose some ecosystem-based adaptation approaches related to water pollution. Finally, in the fourth section, we discuss how the grey water footprint for nitrate or total phosphorous could be an EbA option for improving decision-making and water security in subtropical catchments under change.

2 Material and methods

2.1 The case study area

Two of the most vulnerable areas in the Brazilian southeast are the upper Tietê (drainage area 7390 km²) and Piracicaba–Capivari–Jundiaí (PCJ ;drainage area 14 178 km²) watersheds, particularly due to their high population: 18 million inhabitants in the upper Tietê River basin and 5 million in PCJ (Sao Paulo, 2017; IBGE, 2010). In an attempt to ensure public water supply, the government built the Cantareira system, an inter-basin transfer, in two stages: (a) between 1968 and 1974, at the end of a 35-year period that underwent a severe drought in the Piracicaba watershed; and (b) in 1982, with the inclusion of two additional reservoirs that regularized the increasing rainfall from the mid-1970s until 2005 (Zuffo, 2015). The study area comprises the part of the Cantareira system that drains into the Piracicaba River and which is the headwater of the Piracicaba basin (Fig. 1). This basin is located on the borderline of the state of Minas Gerais and São Paulo. This part of the water supply system, in the Piracicaba watershed, consists of three main reservoirs, named after the rivers, damming the Jaguari–Jacareí, Atibainha and Cachoeira watersheds (drainage areas are 1230, 392 and 312 km², respectively). These rivers are main tributaries of the Piracicaba River, which is a tributary of the Tietê River system on the left bank of the Paraná Basin. The Cantareira system consists of two more reservoirs out of the Piracicaba River basin, Paiva Castro and Águas Claras, which are not part of our study area.

Figure 1 Location of the Cantareira water supply system in the Piracicaba and upper Tietê watersheds.

With respect to the water quality, the headwaters of the Cantareira system are classified as class 1 (best quality) for the Jacareí, Cachoeira and Atibainha watersheds and class 2 for the Jaguari watershed, according to the CONAMA Resolution no. 357/2005 (Brazil, 2005) and São Paulo Decree no. 8468/1976 (Sao Paulo, 1976), which means that, with the exception of the Jaguari watershed, the others can be used for supply with only a simple treatment. Regarding the water volume, this region has been intensely impacted by a severe and recent drought (Taffarello et al., 2016a; Escobar, 2015; Whately and Lerer, 2015; ANA, 2015; Porto and Porto, 2014). As a result of this serious water crisis, a new water policy on the average flow of the transfer limits of the Piracicaba watershed to the upper Tietê watershed was postponed from 2014 to May 2017 (ANA, 2015). The Cantareira system is located in the Atlantic Forest biome, which is considered a conservation hotspot because of its rich biodiversity. In spite of that, 78 % of the original forest cover of the Cantareira watershed has been deforested over the past 30 years (Zuffo, 2015). In 2014, the native forest cover was 10 % in Extrema, 12 % in Joanópolis and 21 % in Nazaré Paulista (SOS Mata Atlântica/INPE, 2015). To counteract deforestation, some environmental and financial tradeoffs have been developed in the Cantareira headwaters to protect downstream water quality and the regulation of water flows. These are ecosystem-based adaptation initiatives, in which rural landowners receive economic incentives to conserve and/or restore riparian forests and implement soil conservation practices. The first Brazilian EbA approach was the Water Conservator (Conservador das Águas) project, created in 2005 and implemented in Extrema, Minas Gerais (Richards et al., 2015, 2017; Pereira, 2013). The Water Producer/PCJ project was developed from 2009 to 2014 in the Cantareira system region (Guimarães et al., 2015), using EbA scenarios and local actions adopting the concept of payment for hydrological ecosystem services (Pagiola et al., 2013; Padovezi et al., 2013) through public–private partnerships, strengthening EbA in Brazil.

2.2 Databases and model adopted

Figure 2 shows the method developed and applied to assess the regulating hydrologic services through grey WF, as well as the spatial data used in this study. The simulations were enhanced by model parameterization with qualitative and quantitative primary data (Mohor et al., 2015a, b; Taffarello et al., 2016b) from six field campaigns between 2012 and 2014, in partnership with ANA, CPRM, TNC-Brazil, WWF, USP/EESC and municipalities. This can reduce uncertainties of the model, facilitate data interpretation and provide consistent information. We installed three data collection platforms (DCP) in catchments at Posses, Cancã and Moinho, as well as level and pressure sensors (see Table 1 and Fig. 8) in paired sub-basins (i) with high original vegetation cover and (ii) in basins that receive payment for ecosystem services (PES) due to participating in the Water Producer/PCJ project.

Figure 2 Methodological scheme for assessing hydrologic services based on greyWF.

Table 1 Sub-basins delimited in SWAT with drainage areas and geographic locations.

^a Indicates new data collection stations installed for experimental monitoring according to ANA/CPRM standards. ^b Indicates experimental stations for research purposes. Source: Taffarello et al. (2016a)

We obtained and organized secondary data from the region upstream of the Jaguari–Jacareí, Cachoeira and Atibainha reservoirs. We then set up a database originating from several sources: Hidroweb (ANA, 2014), Water Sanitation Company of the State of São Paulo (SABESP), Integrated Center for Agrometeorology Information (CIIAGRO, 2014), Department of Water and Power (DAEE), National Institute of Meteorology (INMET) and the Center for Weather Forecasts and Climate Studies (CPTEC/INPE). Supplement Table S1 summarizes all hydrologic, pedological, meteorological and land-use data used as input for the delineation and characterization of the watersheds. The topographical data used was the digital elevation model ASTER Global DEM, second version, 30 m (Tachikawa et al., 2011), available free of charge at http://gdex.cr.usgs.gov/gdex/ (last access: 9 August 2015). The changes in hydrologic services can be evaluated by a large number of models (Carvalho-Santos et al., 2016; Duku et al., 2015; Quilbé and Rousseau, 2007), especially those more user-friendly for stakeholders and policy makers. Simulations in this watershed-scale ecohydrologic model (Borah and Bera, 2004) allow for the quantification of important variables for ecosystem services analysis and decision-making. Some examples of ecohydrologic models with progressive applications in Brazilian basins are the Soil and Water Assessment Tool (SWAT; Bremer et al., 2016; Francesconi et al., 2016; Bressiani et al., 2015), the models reviewed by de Mello et al. (2016), Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) (Sharp, 2016) and Resource Investment Optimization System (RIOS) (Vogl et al., 2016).

The Texas A&M University Soil and Water Assessment Tool (SWAT-TAMU; Arnold et al., 1998; Arnold and Fohrer, 2005) is a public domain conceptual spatially semi-distributed model, widely used in ecohydrologic and/or agricultural studies at a river basin scale (Krysanova and Whyte, 2015; Krysanova and Arnold, 2008). It divides the basin into sub-basins based on an elevation map and the sub-basins are further subdivided into hydrologic response units (HRUs). Each HRU represents a specific combination of land use, soil type and slope class within the sub-basin. The model includes climatic, hydrologic, soil, sediments and vegetation components, transport of nutrients, pesticides, bacteria, pathogens, BMP and climate change in a river basin scale (Raghavan Srinivasan, personal communication, 2014; GASSMAN et al., 2014; Arnold et al., 2012). There have been at least 2600 published SWAT studies (SWAT Literature Database, mid-2016). In the SWAT Purdue conference, held in 2015, 118 studies were presented, of which only 8 % assessed the transport of nutrients in watersheds (SWAT Purdue, Book of Abstracts, 2015). Research using SWAT was developed not only for quantity but also for water quality and ecosystem service assessments (Francesconi et al., 2016; Abbaspour et al., 2015; Duku et al., 2015; Gassman et al., 2014; Fukunaga et al., 2015) and also as an educational tool for comparing hydrologic processes (Rajib et al., 2016).

2.3 Model setup

The initial model setup used the ArcSWAT interface, integrated to ArcGIS 10.0 (Environmental Systems Research Institute). Discretization in sub-basins was carried out, where possible, at the same NCE sites of field investigations. The delimitation of the basin using ArcSWAT requires a drainage area threshold, determined to be 7.1 km², dividing the geographical space to represent the 17 sampling sites in the research field as sub-basins, plus the limits of the three reservoirs'drainage areas, which resulted in 20 sub-basins (Table 1 and Fig. 1b). We highlight that the basin was designed up to the confluence of the Jaguari and Atibaia rivers, forming the Piracicaba River, to integrate all areas of interest in the same SWAT project. The definition of the HRU was carried out using soil maps of the state of São Paulo. (Oliveira, 1999) and land-use maps were developed by Molin (2014) and Molin et al. (2015) from Landsat 5 Thematic Mapper imagery for 2010, using a 1:60 000 scale. The procedure defined 49 HRUs inside the 20 sub-basins, i.e., 49 different combinations of soil type, soil cover and slope classes in our study area. Next, we adapted the land-use map developed by Guimarães et al. (2015), which represents a 2010 land-use scenario for the Cantareira system restoring the most fragile degraded parcels (greatest potential for sediment production), to agree with the land-use classes of Molin (2014). Additionally, we assumed that the second scenario of Guimarães et al. (2015), who used the INVEST model to provide the ecological restoration benefits in the Cantareira system, could be achieved in year 2035, considering the investments provided in the PCJ river plan (Cobrape, 2011) to recover riparian forests. It is worth mentioning that, in the PCJ basin plan, it is called the “trend scenario”. Since in the region the restoration of riparian forests is mostly due to Water-PES projects, which was recognized as an ecosystem-based adaptation (according to CBD, 2010; BFN/GIZ, 2013; Taffarello et al., 2017), we identify the third scenario as S2 + EbA. Thus, Fig. 3 shows the land-use changes over time. In the trend scenario (PCJ-COBRAPE, 2011), the municipalities covered by the Cantareira system could reach a 98 % collection rate with a sewage treatment rate of 100 % and BOD_5,20 removal efficiency of 95 % (PCJ-COBRAPE, 2011). Some studies have suggested including other parameters such as dissolved oxygen, nitrate and phosphate polluting loads, as well as sediments to assess the water quality (Cunha et al., 2014). Regarding the treatment costs for drinking water supply, ecosystem-based adaptation options, such as watershed restoration, seem to be more cost-effective than many technologies for water treatment (Cunha et al., 2016).

Figure 3 Land-use change during 1990 (scenario S1), 2010 (scenario S2) and 2035 (scenario S2 + EbA) in the headwaters of the Cantareira water supply system.

2.4 Calibration and validation

We used the SWAT CUP 5.1.6.2 interfaces and sequential uncertainty fitting (SUFI-2) algorithm for calibrating the quantity and quality parameters and also for validating the simulations in the sub-basins. Quantitative calibration was performed in stations that had more than 2 full years of observed data, i.e., eight stations, namely Posses outlet, F23, F24, F25B, F28, Atibainha Reservoir, Cachoeira Reservoir, Jaguari and Jacareí Reservoir (Table 2). A common test period for all LULC scenarios was selected; in our case, the test period ranged from 1 January 2006 to 30 June 2014. This period has the rain anomaly of drought conditions from 2013 to 2014. The calibration period was from October 2007 to September 2009, the only period with observed data in all of the above eight stations. Validation took place from January 2006 to September 2007 and from October 2009 to June 2014. Calibration and validation of SWAT at the stations with over 2 years of data were rated as good, according to the classification by Moriasi et al. (2007), since the Nash–Sutcliffe efficiency (NSE) criterion (Nash and Sutcliffe, 1970) was greater than 0.65, except for the Posses outlet, which presented the logarithmic Nash–Sutcliffe efficiency of the logarithm of discharges (NSElog) (using the logarithm of streamflow, a criterion that gives greater weight to smaller flow rates) of less than 0.5, rated as unsatisfactory. The percent bias (Pbias) statistic indicates the bias percentage of simulated flows relative to the observed flows (Gupta et al., 1999). Thus, when the Pbias value is closer to zero, it results in a better representation of the basin and in lower estimate tendencies (Moriasi et al., 2007). As a general rule, if $\left|\text{Pbias}\right|<\mathrm{10}$ %, it means a very good fit; 10 % $<\left|\text{Pbias}\right|<\mathrm{15}$ % means good; 15 % $<\left|\text{Pbias}\right|<\mathrm{25}$ % means satisfactory; and $\left|\text{Pbias}\right|>\mathrm{25}$ % means the model is inappropriate. On the other hand, the NSE coefficient translates the application efficiency of the model into more accurate predictions of flood flows, using the following classifications: NSE >0.65 the model is rated as very good; 0.54 < NSE <0.65 the model is rated as good; and between 0.5 and 0.54 it is rated as satisfactory.

Table 2 Characteristics of quantitative calibration and validation of SWAT in studied catchments (Moriasi et al., 2007). Area delimited by the digital terrain model (adapted from Mohor and Mendiondo, 2017).

In the results obtained for different basin scales (Fig. 4), the Pbias and NSE coefficients (including NSE of logarithms) indicate adequate quantitative adjustments. As the SWAT simulations include more than 200 parameters, based on research from the literature (Duku et al., 2015; Bressiani et al., 2015; Arnold et al., 2012; Garbossa et al., 2011), we selected approximately 10 parameters (see Table 3) to complete the calibration to simulate streamflow processes and nutrient dynamics. These parameters refer to key processes which represent soil water storage, infiltration, evapotranspiration, flow channel, boundary conditions (see Mohor et al., 2015b) and main water quality processes at hillslopes. Although our calibration is mainly focused on water yield as total runoff, freshwater quality features through pollutant loads were performed in the scenarios. Further comments related to the existing literature for selected model parameters are depicted in Supplement Sect. S3 with comments on sensitivity analysis to select model parameters used in this paper (the Supplement).

Figure 4 Model calibration related to drainage areas of catchments in the Cantareira system.

Table 3 Calibrated SWAT parameters in the headwaters of the Cantareira water supply system. SCS is related to the method developed by the Natural Resources Conservation Service of USDA, previously named the Soil Conservation Service. USLE is the universal soil loss equation.

Moreover, to reduce the uncertainty of our predictions, we used approximately 2500 primary data derived from an earlier stage of this research (Taffarello et al., 2016a). As a parameterization result of field investigations and ecohydrologic modeling, Fig. 5 shows parts of the calibrated model performance (lines) against field observations (dots with experimental uncertainty) for flow discharges, nitrate and total phosphorus loads for catchment areas ranging from 7.1 to 508 km². Finally, other water quality variables were studied based on data from field sampling.

Figure 5 Comparison between flow discharges (upper part) and nitrate loads (lower part), through observed flows (dotted lines), simulated by SWAT (solid lines) and field validation through instantaneous experimental samples (marked points with uncertainty intervals) at monitored stations of Posses outlet (left part), F23 Camanducaia (center part) and F24 Cachoeira (right part). The uncertainty bars were determined using instantaneous velocities measured in the river cross-sections during 2013/14 field campaigns (Taffarello et al., 2016a). The uncertainty bars represent the minimum and maximum values of measured streamflow and pollutant loads in a cross section of the river during a field campaign of headwater catchments.

We highlight some SWAT model limitations when we compare the simulated to observed water flows, especially in the dry season. For example, when the model was discretized on a daily resolution, the adherence level between the observed and simulated flows was considered good. However, the model did not fit well to the observed values only during a specific time interval of the drought period (i.e., February 2014 to May 2014). These differences were more significant for water quality parameters, such as nitrate and total phosphorous. We point out that the macronutrient loads found in May 2014 were clearly higher than the loads we found in previous sampling, which occurred in wetter periods (Taffarello et al., 2016). For the sample collected in May 2014, the model significantly underestimated the pollutant loads of nitrate. This behavior, arising from the recent and most severe drought faced by the Cantareira system (Nobre et al., 2016; Marengo et al., 2015; Taffarello et al., 2016a; Escobar, 2015; The Economist, 2015; Porto and Porto, 2014), shows the need for the improvement of the SWAT model performance, especially to capture nonlinearities having impacts on regulating ecosystem services during extreme flows. For EbA scenarios, we planned to set up field investigations and SWAT calibrations (see Fig. 5) using the extreme conditions of the 2013–2014 drought through freshwater quality monitoring at the headwaters of the Cantareira system (see Tafarello et al., 2016a).

2.5 The scenarios and a new index for hydrologic service assessment

Differences in flow rates and water quality (for the variables of nitrate, phosphate, BOD_5,20, turbidity and faecal coliforms) for the 20 sub-basins were evaluated using flow-and-load duration curves for the three scenarios proposed in this study: (i) recent past scenario (S1), including the recorded past events for land use in 1990; (ii) current land-use scenario (S2), which considered land uses for the 2010–2015 period as the baseline; and (iii) future land-use scenario (S2 + EbA), supposing a forest cover conversion in the protected areas, through EbA options, according to the PCJ river basin plan by 2035. Using these curves, from the methodology shown by Hoekstra et al. (2011), and based on Duku et al. (2015) and Cunha et al. (2012), we estimated the grey water footprint. Next, we developed a new ecohydrologic index to assess the regulating hydrologic services in relation to the greyWF. This new indicator encompasses the former theory related to environmental sustainability of the greyWF, according to Hoekstra et al. (2011). In this study, as a relevant local impact indicator, Hoekstra et al. (2011) proposed to calculate the water pollution level (WPL) within the catchment, which measures the degree of pollution. WPL is defined as a fraction of the waste assimilation capacity consumed and calculated by taking the ratio of the total of greyWF in a catchment (∑WFgrey) to the actual runoff from that catchment (R_act), or, in a proxy manner, the water yield or mean water yield or long-term period (Q_lp). This assumption is that a water pollution level of 100 % means that the waste assimilation capacity has been fully consumed. Furthermore, this approach assumes that when WPL exceeds 100 %, environmental standards are violated, such as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\text{WPL}={\displaystyle \frac{\sum \text{WFgrey}\left[x,t\right]}{R_{\text{act}}[x,t]}}.\end{array}$$

It is worth mentioning that for some experts, the aforementioned equation can overestimate the flow necessary to dilute pollutants. For that reason, new insights of composite indicators or thresholds are recommended, as follows. The above assumption could overestimate WPL because it would fail to consider the combined capacity of water to assimilate multiple pollutants (Hoekstra et al., 2012; Smakhtin et al., 2005). Conversely, in this study, we define an alternative indicator related to the three following fundamentals. First, the WPL should be extended to a composite index, thereby representing weights of each pollutant related to the actual runoff, here as a proxy of long-term runoff, i.e.,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{WPL}}_{\text{composite}}\left[x,t\right]={\displaystyle \frac{\sum \left\{w\left[x,t\right]\cdot \text{WFgrey}[x,t]\right\}}{R_{\text{act}}[x,t]\cong Q_{\text{lp}}[x,t]}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }\sum w\left[x,t\right]=\mathrm{1},\mathrm{0}\le w[x,t]\le \mathrm{1}.\end{array}$$

For this new equation, weights should be assessed, either from field experiments or even from simulation outputs. Second, we define a threshold value of WPL_{composite,ref} regarding the reference catchments in non-developed conditions which suggest more conservation conditions among other catchments of the same region, as WPL_reference. For this study, we selected the Domithildes catchment as the reference catchment with conservancy measures. From this reference catchment, we define the composite reference index for the water pollution level as WPL_{composite,ref} and, derived from it, the hydrologic service index (HSI), as a non-dimensional factor of comparison between WPL for reference and non-reference catchments, as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\text{HSI}\left[x,t\right]\text{greyWF}={\displaystyle \frac{\text{WPL}\left[x,t\right]-{\text{WPL}}_{\text{composite,ref}}}{{\text{WPL}}_{\text{composite,ref}}}}.\end{array}$$

3 Results

In the following section, we present the results from field observations, which are useful not only for ecohydrologic parameterization, but also for elucidating features regarding greyWF and hydrologic services. Next, we compare the water yield and greyWF outputs from simulations under LULC scenarios, including EbA options, to finally propose a new hydrologic services indicator.

3.1 Data from field sampling

Some of the water quality and quantity variables from our freshwater monitoring are useful to assess the hydrologic services; thus, they are presented in Table 4. These variables were selected due to their relationship with anthropic impacts on the water bodies and because of their importance for sanitation. Among the water quality variables sampled in the field step of the research (see Taffarello et al., 2016a, b), we highlight turbidity because it indicates a proxy estimation about the total suspended solids in lotic environments (UNEP, 2008), related to the LULC conversion, and reflects the changes in the hydrologic services. Figure 6 shows the direct correlation between turbidity and size of the sub-basins. Turbidity can indirectly indicate anthropic impacts in streams and rivers (Martinelli et al., 1999). The lower turbidity mean values were observed in two more conserved sub-basins (which presented higher amounts of forest remnants): 2 NTU in the reference Cancã catchment (Domithildes) and 5 NTU in the upper Posses. Otherwise, we found a positive relationship between nitrate concentrations and both discharge and mean water level (Fig. 7). It can be inferred that higher concentrations of macronutrients would be found in downstream areas. This trend can be associated with the nutrient migration (Cunha et al., 2013) and land-use change (Zaffani et al., 2015), as well as point source pollution. In addition, the absence of the riparian forest in 70 % of the protected area (36.844 ha) of the Cantareira system (Guimarães et al., 2015) can increase the sediment transport from riparian areas to rivers and make pollutant filtration more difficult, leading to higher nitrate concentrations downstream.

Figure 6 Experimental sampling of turbidity (size of circles), observed flows and mean velocities in river cross sections of 17 catchments in Cantareira system headwater (October 2013–May 2014). This illustration shows the high interdependence and complexity in integrating any standard parameterization, at a regional scale, of the SWAT model, linking potential scenarios of LULC, water yield and freshwater quality in medium-size basins and headwaters.

Figure 7 Multidimensional chart of hydraulic and water quality variables sampled in field campaigns in the headwaters of the Cantareira water supply system between October 2013 and May 2014.

Figure 8 Study area divided into sub-basins for hypothesis testing using semi-distributed SWAT model.

3.2 LULC change scenarios

The variations in LULC affect freshwater quality, which, in turn, affect the dynamics of aquatic ecosystems (Zaffani et al., 2015; Hamel et al., 2013; Bach and Ostrowski, 2013). These changes impact the hydrologic services, especially regulating and supporting ecosystem services (Mulder et al., 2015; Molin et al., 2017). The LULC of each sub-basin, according to a past-condition scenario (S1, in 1990), a present-condition (S2, in 2010) and a future (S2 + Eba, in 2035) LULC scenario, using the same weather input data files, is shown in Table 5. We evaluated the effects of LULC change scenarios in 20 catchments in the Jaguari, Cachoeira and Moinho sub-basins, southeast Brazil (Fig. 8). Concerning the land-use change, the main soil uses 25 years ago were pasture (in 50 % of the sub-basins) and native vegetation (in 45 % of the sub-basins). According to Molin (2014), the 5 % of the remaining area was divided into vegetables, eucalyptus, sparse human settlements, bare soil and mining. The main activity in the past (1990) was extensive cattle raising for milk production by small producers in the region (ANA, 2012; Veiga Neto, 2008). By assessing the temporal trends of increment or reduction of native remnants, we examined the periods 1990–2010 vs. 2010–2035. From 1990 to 2010, the percentage of forest increased by 50 % in the Domithildes sub-basin, which was the reference catchment of the Water Producer/PCJ project (see Taffarello et al., 2016a), Moinho, Cachoeira dos Pretos, F34, B. Jacareí, B. Atibainha, B. Cachoeira, Pq Eventos, F25B and B. Jaguari (Fig. 9). Concerning the period from 2010 to 2035, the model was set up considering an increase in native vegetation in all sub-basins from forest remnants in 2010 and from the new BMP of reforestation with native species in 20 sub-basins by 2035 (Fig. 9). The hydro-services in the Posses and Salto catchments and in the Cachoeira sub-basin will be increased by 2035 as a function of the efforts on EbA which currently exist in the region (Richards et al., 2015, 2017; Santos, 2014).

Figure 9 Native forest cover in S1 (1990), S2 (2010) and S2 + EbA (2035).

3.3 Water yield as a function of soil cover

In this research, we chose to use qualitative and quantitative duration curves for the integrated assessment of availability and quality of water. The flow-and-load duration curve, comparable to histograms of relative cumulative frequencies of flows and loads of a water body, is a simple and important analysis in hydrology (Collischonn and Dornelles, 2013). In quantitative terms, the flow duration curve shows the probabilistic temporal distribution of water availability (Cruz and Silveira, 2007), relating the flow in the river cross section to the percentage of time in which it is equalled or exceeded (Cruz and Tucci, 2008).The three scenarios S1, S2 and S2 + EbA resulted in different flow values for the 20 sub-basins (Fig. 10). Based on the arithmetic mean of time series of monthly water yields, related to catchment areas, and assessed for all modeled sub-basins (N=20), the results show average values of water yield: 31.4±25.2 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ for S1 (1990), 14.9±11.5 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ for S2 (2010) and 21.4±15.3 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ for S2 + EbA (2035). This very high variation can be due to the complexity of river basin systems and the various sources of uncertainty in the representation of ecohydrologic processes.

Figure 10 Flow duration curves under three LULC scenarios: S1(1990), S2(2010) and S2 + EbA(2035) at headwaters of the Cantareira water supply system.

The three analyzed scenarios and the ecohydrologic monitoring provide different types of information for the same catchments. The 52 % decrease in water yield between S1 (1990) and S2 (2010) scenarios, as ($(\mathrm{14.9}-\mathrm{31.3})/\mathrm{31.3}\times \mathrm{100}$), might be related to a marginal increase in the eucalyptus cover. In fact, from 1990 to 2010, eucalyptus cover increased +6.8 % in total land cover, but +181 % in relative terms. Another possible explanation is the decrease in native vegetation from 1990 to 2010, with −1.8 % in total land cover, but −4.3 % in relative terms. In parallel, we evaluated the water yield. Thus, the flow-and-load duration curves summarize the flow and pollutant load variability, thereby showing potential links and impacts for aquatic ecosystem sustainability (Cunha et al., 2012; Cruz and Tucci, 2008). From these curves, we obtained two different behaviors for the studied sub-basins (Fig. 10).

Behavior I: the water yield in 2010 was reduced in relation to 1990 and the water yield in 2035 might exceed the 1990 levels. The examples are as follows: the Upper Jaguari (Alto Jaguari), Cachoeira sub-basin (including the Cachoeira dos Pretos, Chalé Ponto Verde, Ponte Cachoeira, F24 outlet) and Moinho catchments.

Behavior II: the water yield after 2010 was reduced until 2035 and this water yield recuperation was not possible for the values in 1990. Examples, in decreasing size of drainage areas, are as follows: Atibainha, B. Jaguari, F25B, Parque de Eventos, F23, B.Atibainha, F34, F30, Salto, Posses outlet, Domithildes, Portal das Estrelas (middle Posses).

On the one hand, according to Fig. 11, the water yield of S1 is inversely proportional to the land use of mixed forest cover. The water yield in S2 indicates a constant value of approximately 17 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$. Moreover, for the S2 + EbA scenario, which incorporates the EbA approach through BMP, the water yield is approximately 17 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$, but with a slight increase in the water yield when the percentage of forest cover is higher than 50 %. Presumably, this slight increase in the water yield would be related to the type of best management practices of the recovery forests, which still did not achieve evapotranspiration rates of the climax stage. In the riparian forest recovery, evapotranspiration rates are lower and, thus, a greater amount of precipitation reaches the soil and rivers through the canopy. This process could benefit other hydrologic components, such as runoff, increasing water flows into the rivers. This effect can possibly explain the behavior I catchments (see Fig. 10).

Figure 11 LULC scenarios for specific water yield for 20 drainage areas at Jaguari, Cachoeira and Atibainha watersheds, according to S1 (1990), S2 (2010) and S2 + EbA (2035) scenarios.

3.4 Relationships between land-use and land-cover change and grey water footprint

For an integrated assessment of hydro-services, we analyzed the spatiotemporal conditions of load production at the sub-basin scale (see more information on Sect. S4 “Comments on differences in land use and land cover in sub-basins studied”, in the Supplement). As we studied rural sub-basins, water pollution was mainly produced by diffuse sources, such as fertilizers and agrochemicals. In this context, we evaluated the evolution of greyWF to show nitrate (N-NO₃), total phosphorus and sediment (Sed) yields (indicated by turbidity) of scenarios S1, S2 and S2 + EbA. First, we calculated the nitrate loads generated from the 20 sub-basins in the three scenarios. Second, we did the same for total phosphorous loads and sediment yields. Third, considering the river regime, we calculated the greyWF for nitrate, total phosphorous and sediments in each sub-basin to develop a new composite index that assesses the sustainability of hydrologic services.

Concerning nitrate, the sampled concentrations were low. In addition, SWAT simulations also produced very low outputs, and the greyWF-NO₃ varied from 0.11 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ (in Atibainha sub-basin in the S2, 2010, scenario) to 2.83 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ (in the middle Posses catchment, Portal das Estrelas, under the S2 + EbA, 2035, scenario). Considering Brazilian water quality standards for nitrate, the maximum allowed concentration is 10 mg L⁻¹ (Brasil, 2005). These low amounts of nitrate loads make the greyWF-NO₃ fall to low values in the three analyzed scenarios (between 1 and 10 %; Fig. 12a). In relation to total phosphorous, the load duration curves from S1, S2 and S2 + EbA scenarios showed disparities. For example, the greyWF-TP decreased in all sub-basins between 1990, 2010 and 2035. From 2010 to 2035, the model predicts a new behavior for the greyWF-TP.

Figure 12 Fraction of water yield (mean Q) compromised by the grey water footprint of nitrate (greyWF-NO3), total phosphorous (greyWF-TP) and sediments (greyWF-Sed) vs. drainage area (a) and vs. selected sub-basins (b).

Results of the greyWF for TP, NO₃ and sediments enabled us to infer hydrological regionalization for nutrient loads. Among the 20 sub-basins studied, we selected 2 sub-basins as study cases to illustrate the links between LULC and greyWF: (1) the Upper Jaguari and (2) Domithildes. The reasons for selecting the two sub-basins among the 20 sub-catchments are detailed in Sect. S5.

3.4.1 Case study I: Upper Jaguari sub-basin

The Upper Jaguari (Fig. 13) has 302 km² and is the second most upstream sub-basin within the Cantareira system (downstream of only the F28 sub-basin, with 277 km²). Comparing scenario 1990 (S1) and 2010 (S2), the results showed evidence that the native forest decayed approx. 10 %. Indeed, scenario 2035 (S2 + EbA) still assumes a very small decrease in the native forest. This decrease may be due to the increase in secondary forests by BMP, which could stabilize the native forest LULC by 70 % until 2035. The mean annual simulated water yields, in spite of high variability of simulated scenarios, pointed out values of 18 L s⁻¹ km² (1990, S1), 13 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ (2010, S2) and 21 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ (for 2035, S2 + EbA).

Figure 13 Synthesis chart of case study of the Upper Jaguari sub-basin (drainage area =302 km²). (a) Localization at the drainage areas of the Cantareira system; (b) LULC conditions for scenarios S1 (1990), S2 (2010) and S2 + EbA (2035); (c) comparison of water yields simulated for conditions of S1, S2 and S2 + EbA; (d) water yield scenarios compared with intra-annual regime of S2 + EbA scenario; (e) comparison of duration curves of flows for S1, S2 and S2 + EbA conditions; (f) duration curves of N-NO3 loads for S1, S2 and S2 + EbA.

3.4.2 Case study II: Domithildes headwater

The Domithildes catchment (9.9 km²) is located in the Cancã catchment. Similar to the Upper Jaguari, Domithildes is one of the most conserved sub-basins, mainly with native forests. The native forest fraction remained constant (see Fig. 14) from S1 (51 % in 1990) to S2 (52 % in 2010). However, unlike the Upper Jaguari sub-basin (see Fig. 13), native vegetation could increase by 56 % in S2 + EbA (2035). Due to the fact that Domithildes was adopted as a reference basin for Water Producer/PCJ, the augmented fraction of native forest by 2035 could show an increase in secondary forest. Regarding water yield, the Domithildes catchment was classified as a second type of sub-basin behavior (Sect. 3.3). There is a positive increment of water yield between 2010 (∼18 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$) and 2035 (∼23 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$), although this situation may not achieve values obtained for S1 conditions in 1990 (∼29 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$).

Figure 14 Synthesis chart of case study the Domithildes catchment (drainage area =9.9 km²). (a) Localization at the drainage areas of the Cantareira system; (b) LULC conditions for scenarios S1 (1990), S2 (2010) and S2 + EbA (2035); (c) comparison of water yields simulated for conditions of S1, S2 and S2 + EbA; (d) water yield scenarios compared with intra-annual regime of S2 + EbA scenario; (e) comparison of duration curves of flows for S1, S2 and S2 + EbA conditions; (f) duration curves of N-NO3 loads for S1, S2 and S2 + EbA.

3.5 Results of a new index for hydrologic service assessment

The new index for hydrologic service assessment was developed as a simple relation between greyWF and water yield, using a fraction between water demand (numerator) and availability (denominator). Some authors commonly use this fraction as a direct approach to water scarcity (i.e., Smakhtin et al., 2005; Hoekstra, 2013; McNulty et al., 2010; among others). Therefore, we first assessed greyWF by respective drainage basins (Fig. 15). Then, we calculated the water pollution levels. The results in Fig. 16 show the composite water pollution level (WPL_{composite,ref}) vs. drainage areas and compared with the HSI. The baseline WPL_{composite,ref} is related to the Domithildes catchment (horizontal, dotted line in Fig. 16). This line divides the graph into two regions: less sustainable basins (HSI>0) and more sustainable basins (HSI≤0). The more sustainable basins (HSI<0) are the Salto; Cachoeira nested catchments (Cachoeira dos Pretos, Chalé Ponto Verde and Ponte Cachoeira); and the F28, F24 and Upper Jaguari basin.

Figure 15 Relationships between grey water footprint for nitrate (a) and total phosphorous (b) according to three LULC scenarios (1990, 2010 and 2035) and size of the drainage areas of headwaters in the Cantareira water supply system.

Figure 16 Hydrologic service index (circle ratio) related to drainage area of river basin (horizontal axis) and composite of water pollution index (vertical axis) for the S2 + EbA scenario: equal weights of nitrate, total phosphorus and dissolved sediments are expressed in WPL_{composite,ref}.

3.6 Comparison of field investigation and modeled scenarios

Field and experimental data (Taffarello et al., 2016a) with modeled scenarios of land-use and land-cover change, including the EbA hypothesis, were integrated into a summary figure in the Supplement (see Fig. S1).

4 Discussion

This section discusses field data, LULC change scenarios, greyWF and water yield, not only in general aspects, but also in selected catchments, which is mentioned in Sect. 3.

4.1 On field data

Other conserved sub-basins also presented low mean values of turbidity (<6.5 NTU), e.g., intervention Cancã catchment (5 NTU) and Cachoeira dos Pretos (6 NTU). We found the highest turbidity above 40 NTU, which is considered the maximum established water quality standard for Brazilian class 1 (BRASIL, 2005), at Parque de Eventos (283 NTU), at F23 (180 NTU) and at Salto outlet (160 NTU). However, these three sampling sites are located at water bodies of class 2, where the maximum turbidity allowed is up to 100 NTU (BRAZIL, 2005). Due to these areas having the highest urbanization among the sampled sites, they are in noncompliance with Brazilian environmental standards. Arroio Jr. et al. (2013) found a decreasing relation between turbidity and drainage areas in another catchment located in São Paulo state. Temporal turbidity patterns show that on the one hand, in 11 out of 17 monitored sites, the higher values of turbidity occurred in December 2013, the only field campaign with significant precipitation (35.3 mm) and with a higher antecedent precipitation index (API =123.7 mm). This can be due to carrying allochthon particles, which are drained into rivers by precipitation. Analogously, Arroio Jr. et al. (2013) also observed higher turbidity in the rainy season (December 2012) which can lead to erosive processes. On the other hand, Zaffani et al. (2015) showed that turbidity did not vary over the hydrologic year in medium-sized, rural and peri-urban watersheds ranging from 1 to 242 km². In this case, other factors may have had an influence, such as deforestation, seasonal variability, soil use type, sewage and mining (CETESB, 2016; Tundisi, 2014).

4.2 On LULC change scenarios

The main soil uses are pasture (58 % of the drainage area of the sub-basins) and forest (40 % of the drainage area of the sub-basins) in 2010. From 1990 to 2010, there was a significant conversion of soil cover, with a small reduction of pasture areas (−2 %) and native remnants (−5 %) and with a progressive increase in eucalyptus, a non-native tree species in Brazil. Eucalyptus soil use varied from an increase of +1 % within Posses up to an increase of +31 % in the Chalé Ponto Verde sub-basin in 2010 (compared to 1990). Eucalyptus cover, however, did not achieve 10 % of the soil uses in any of the simulated sub-basins in 1990. In the third scenario (S2 + EbA), we hypothesized incentives of public policies for forest conservation and restoration, due to the strengthening of EbA in the Cantareira system. This could lead to an increase in native vegetation reaching percentages of 15 % in the Posses outlet and 69 % in the F28 sub-basin. In this scenario, the higher percentages of native vegetation would occur in the sub-basins F28, Upper Jaguari and Cachoeira dos Pretos. Despite this general increase in native forest cover, we highlight the deforestation which occurred in the F23 sub-basin in the Camanducaia River. Currently, although the basin has 34 % of native forest cover, this rate has tended to decrease since 1990. The F23 outlet (sub-basin 2) had 37 % of native forest cover in 1990, which then became 34 % in 2010, and the S2 + EbA scenario predicts that F23 could reach 36.2 % of native forest by 2035, returning to the percentages found in 1990. Another critical situation is the Posses outlet (SWAT sub-basin 6). Despite the conservation efforts which have been made in the region through the Water Conservator project (see Richards et al., 2015; Santos, 2014; Pereira, 2013), the current percentage of native remnants is 13 %, which may be 16 % in 2035, however, not achieving the rate in 1990 (22 %). This can potentially disrupt the regulation of hydrologic services provided by Posses sub-basin and needs to be evaluated in depth. Spatiotemporal patterns of the main soil uses which compete with forest cover are analyzed, e.g., pasture and eucalyptus. First, related to pasture, it can be observed that pasture was the main use in the past in 60 % of the sub-basins (in 1990) and, currently, it has become the main LULC – approximately 40 %. Our scenarios indicate that due to EbA strengthening, encouraging the links between environmental conservation and forest restoration, 20 % of the sub-basins could be mainly occupied by pasture (sub-basins 2, 4, 6 and 7). This rate is reasonable, considering rural sub-basins. Moreover, the reduction in pasture in the Cantareira system was more evident in the 1990–2010 period than in the 2010–2035 scenario. This can be explained by, at least, three factors: (i) rural landowner awareness of the relevance of converting pasture to native forest to generate and maintain ecosystem services in the Cantareira system (Extrema, 2015; Padovezi et al., 2013; Gonçalves, 2013; Veiga-Neto, 2008; Richards et al., 2015, 2017); (ii) seasonal changes in the ecosystem structure which can increase the ecosystem resilience (Mulder et al., 2015); and an observed significant increase, mainly in the 1990–2010 period, of nonnative species plantations. Second, regarding the eucalyptus cover, the future scenario shows an increasing threat to the regulating and supporting services as a result of the exotic forest expansion. In 2035, eucalyptus cover may include, on average, 12 % of the total area of the 20 catchments studied here. This is significant in comparison with the 10 % in 2010 and only 2 % in 1990 for the same catchments. The scenario for 2035 shows that the maintenance of hydrologic services deserves attention, because eucalyptus monoculture can potentially impact not only the headwaters, but entire landscapes, threatening the ecosystem dynamics. Moreover, these plantations, with an average wood yield of 50 to 60 m³ of Eucalyptus urograndis per hectare, need high quantities of agrochemicals, due to the low diversity of the population and low adaptation to climate change (Kageyama and dos Santos, 2015). In short, here we highlight the threat to biodiversity that has been brought by alien species in headwaters and the changes that it can promote in native species (Hulme and Le Roux, 2016) which, in turn, impact the ecosystem services.

4.3 On water yield and LULC

On the other hand, we observed in the Posses, Salto, Jaguari, Cancã and Atibainha catchments an inverse situation (behavior II). This effect can be related to the hydrologic response produced by (a) type of catchment; (b) size of catchment; (c) the low soil moisture in the red-yellow latosol (Santos et al., 2006), which did not favor high evapotranspiration rates; (d) the riparian forest, originating from the EbA or Water-PES actions, that should still be at the initial stages, not achieving a climax in 20 years (this explanation therefore assumes that the baseline of PES actions was in 2015, although there are examples of restored forests in Extrema, Minas Gerais, with high evapotranspiration rates, as can usually be found in climax forests); and (e) unpredictability, nonlinearity and uncertainty (Pahl-Wostl et al., 2007; Lima and Zakia, 2006). Riparian native forests, eucalyptus and riparian forests in recuperation (shown here as orchard) have different hydrologic responses. There is still a lack of knowledge regarding the influence of different types and phases of vegetation on the hydrologic processes. Bayer (2014) found that the vegetation height and leaf area index are inversely proportional to the water flows, which corroborate previous studies (Hibbert, 1967). Riparian forest restoration increases the mean evapotranspiration, reducing the water yield (Molin, 2014; Salemi et al., 2012; Lima and Zakia, 2006; Andréassian, 2004). Restoration increases the water storage capability into the catchment throughout the riparian zone, contributing to the higher water flow in the dry season (Lima and Zakia, 2000). This can lead to unexpected results regarding water yield. Furthermore, at small catchments of temperate climate, researchers estimated that deforestation in 40 % of the catchments would increase the runoff of 130±89 mm year⁻¹ considering the entire water cycle in the catchment scale (Collischonn and Dornelles, 2013). In addition, there is high dispersion in the results-based monitoring (usually, in paired catchments or nested catchment experiment), which makes it more difficult to predict the flow as a result of soil use conversion. Similarly, we found high dispersion in the comparison between water yields vs. different land cover in 20 sub-basins of the subtropical climate (Fig. 11). Moreover, BMPs have been in progress since 2005 in the Posses outlet (sub-basin 6, Table 5) and middle Posses (Portal das Estrelas, No 7) and since 2009 in the Domithildes, F30 and Moinho catchments (sub-basins 9, 11 and 20, respectively). These BMPs originated from the Water Conservator and Water Producer/PCJ projects. In these cases, we recommend that public agencies take care when defending PES as inductors of more water availability (ANA, 2013). Parts of these results and previous investigations, which were made through NCE (Taffarello et al., 2016a), point out the opposite; i.e., in the more conserved catchments, we found lower water yields. Despite the fact that there are many Water-PES programs in Brazil (Pagiola, von Glehn and Taffarello, 2013; Guedes and Seehusen, 2011), measurements of the effect on water yield under forest restoration are still lacking in tropical and subtropical conditions (Taffarello et al., 2016a; Salemi et al., 2012). However, the benefits of riparian forests on water quality, margin stability, reduction of water erosion and silting are clear in the scientific literature (Santos, 2014; dos Santos et al., 2014; Studinski et al., 2012; Udawatta et al., 2010).

4.4 On greyWF, LULC and water yield in selected catchments: Upper Jaguari and Domithildes

The discussion of the variability in greyWF and water yield is based on the hydrologic conditions simulated in the test period from 2006 to 2014. In turn, this test period was selected due to high availability of rainfall stations under operation, which would potentially better perform distributed modeling at several sub-basins using SWAT. For the three scenarios simulated, the relationships between the native forest cover and mean water yield are different from each other. On the one hand, in the Upper Jaguari, for scenario S1 (1990), the higher the native forest cover, the lower the water yield. This scenario behavior is extended at experimental sites and even strongly documented in the literature (Salemi et al., 2012; Collischonn and Dornelles, 2013; Ellison et al., 2012). For scenario S2 (2010) the water yield seems not fully related to native forest LULC, oscillating around an average value of 18 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$. In scenario S2 + EbA (2035), however, there is a slight increase in water yield when native forest cover is higher than 50 %. This proportional relation between water yield and forest cover in the S2 + EbA is both controversial and contrary to results published by some authors (e.g., Collischonn and Dornelles, 2013; Salemi et al., 2012). For example, monitoring data shows a reduction in the water yield with higher native forest land cover (Taffarello et al., 2016a). Salemi and coauthors, in a review on the effect of riparian forest on water yield, found that riparian vegetation cover decreases water yield on a daily to annual basis. Furthermore, the greyWF-NO₃ of the Upper Jaguari basin showed 0.14 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ for scenario S1 (1990), increased to 0.23 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ for scenario S2 (2010) and could grow to ca. 0.54 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ in S2 + EbA scenario (in 2035). However, this result is different from the one expected in the hypothesis testing through modeling. The null hypothesis states that increasing native forest cover is correlated to decreasing nutrient loads flowing to streams. The results, modeled by SWAT, predicted an increase in the greyWF by 2035. The simulated increase in the native forest (approx. +5 %) appears to be insufficient for buffering nitrogen loads from animal excrements such as mammals or zooplankton. For a more in-depth analysis, other factors that influence the greyWF should be evaluated thoroughly.

Additionally, in the Domithildes catchment (reference catchment), other factors, such as native vegetation, could influence the hydrologic cycle decreasing water yields in the 2010 scenario (S2). One explanation of this water yield decrease could be the positive LULC of Eucalyptus sp. to +5 % in 2010 (S2). Regardless of other factors, and increase of 1 % of the eucalyptus land-use fraction in Domithildes will represent a decrease of 2 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$ of water yield, or −63 mm per year, in the same range of results reported by Salemi et al. (2012) and close to Smethurst et al. (2015).

Comparing seasonal water yields, the results showed higher variability around monthly flow averages for the S2 + EbA (2035) scenario. These deviations in monthly flows by the S2 + Eba (2035) scenario were higher in wetter months between November and March. The regulation of water yield, in both rainy and dry conditions, is more effective when quantified through variance (Molin, 2014). In spite of these uncertainties, scenarios modeled by SWAT estimated the highest mean monthly water yield in February (38 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$) and the lowest mean monthly water yield in September and October (8 $\mathrm{L}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{2}}$). On the one hand, the results showed that a growing rate of native vegetation LULC since 2010 would serve to attenuate both environmental flow peaks, especially in the rainy season (see flow duration curves), and pollutant filtration (see duration curves of N-NO₃ loads).

On the other hand, the more native forest cover, the lower the water yield (Bayer, 2014; Molin, 2014; Burt and Swank, 1992). Thus, the progressive increase in water yield from 2010 to 2035, compared to a higher total forest cover, could indicate other factors, such as forest connectivity, forest climax and secondary factors, such as BMP, that could produce nonlinear conditions of water yield from the local scale to the catchment scale.

5 Conclusions and recommendations

Although the water–forest system interaction is a classic issue in hydrology, the impacts of vegetation on qualitative and quantitative aspects of water resources need to be better understood. Supported by field experiments and qualitative and quantitative simulations under different scenarios including EbA options with BMP, our results showed evidence of nonlinear relationships among LULC, water yield, greyWF of nitrate, total phosphorus and sediments, which irreversibly affect the composite of the water pollution level, the definition of WPL of reference (here established at the Domithildes catchment) and the hydrologic service index. Although there was a coherent and proportional relation between the observed mean river velocity and observed specific flow, experimental evidence still depicted outliers, not only in reference catchments with EbA/Water-PES options, but also in intervention catchments with no EbA/Water-PES options. This evidence points out illustrative examples of how complex LULC options from EbA would be exhaustively sensed into hydrological parameters and simulation scenarios using SWAT or other distributed models. Despite using a semi-distributed model for assessing non-point sources of pollution mainly tested under different LULC scenarios, our results showed that the intrinsic nature of flow–load duration curves, LULC and greyWF are constrained to high uncertainties and nonlinearities both from in situ sampling and from processes' interactions of modeling. Our results show the need to evaluate many uncertainty sources, such as model sensitivity analysis, observed streamflow data, ecohydrologic model performance and residual analysis. To attain the goals of EbA, using HSI through greyWF assessment and the composite of WPL, some conditions are needed to better fit models to field observations, as follows: (i) monitoring and, if possible, constraining illegal inputs of highly concentrated pollutants, especially from growing urban settlements; (ii) restoring riparian vegetation, especially at HRUs where EbA scenarios introduce more sensitivity of water yields and greyWF; and (iii) modeling EbA effects at HRUs where trapping and removing inflowing sediments are more evident. For the assessment of the health of river ecosystems, we used HSI, flow regimes and the WPL_{composite,ref}, which could bring an alternative method of environmental flow recommendations. Although the role of vegetation on streamflow has been widely studied, very few investigations have been reported in Brazil with control nutrient sources, transportation and delivery. Moreover, further field and modeling research is needed when integrating LULC, EbA and greyWF through hydrologically distributed models. Thus, future research could clarify the influence of vegetation on water quality and the role of anthropogenic and natural drivers in ecohydrologic processes on a catchment scale.