Multidecadal change in streamflow associated with anthropogenic disturbances in the tropical Andes

Andean headwater catchments are an important source of freshwater for downstream water users. However, few long-term studies exist on the relative importance of climate change and direct anthropogenic perturbations on flow regimes in these catchments. In this paper, we assess change in streamflow based on long time series of hydrometeorological data (1974–2008) and land cover reconstructions (1963– 2009) in the Pangor catchment (282 km) located in the tropical Andes. Three main land cover change trajectories can be distinguished during the period 1963–2009: (1) expansion of agricultural land by an area equal to 14% of the catchment area (or 39 km) in 46 years’ time, (2) deforestation of native forests by 11 % (or−31 km) corresponding to a mean rate of 67 ha yr, and (3) afforestation with exotic species in recent years by about 5 % (or 15 km). Over the time period 1963– 2009, about 50 % of the 64 km of native forests was cleared and converted to agricultural land. Given the strong temporal variability of precipitation and streamflow data related to El Niño–Southern Oscillation, we use empirical mode decomposition techniques to detrend the time series. The long-term increasing trend in rainfall is remarkably different from the observed changes in streamflow, which exhibit a decreasing trend. Hence, observed changes in streamflow are not the result of long-term change in precipitation but very likely result from anthropogenic disturbances associated with land cover change.


Introduction
Andean headwater catchments are an important source of freshwater for downstream water users (Urrutia and Vuille, 2009;Roa-García et al., 2011). Although the ecosystems in the tropical Andes have been modified by anthropogenic disturbances for at least 7000 years (Bruhns, 1994), it is only since the early 20th century that natural habitats have undergone extensive transformation (White and Maldonado, 1991). The demand for agricultural land has led to an expansion of the agricultural frontier at the expense of natural ecosystems. The magnitude and intensity of land use change has increased rapidly from the second half of the 20th century, as a result of population growth, socio-economic development, rural-urban and international migration, and land reform programs (Vanacker et al., 2003;Grau and Aide, 2007;Curatola Fernández et al., 2015). The agrarian land reforms of the 1960s and 1970s led to a redistribution of the land ownership but also promoted rapid colonization of so-called vacant lands, which were often covered by native forests (Balthazar et al., 2015).A major concern for sustainable development in the tropical Andes is the increasing demand for freshwater ecosystem services (Harden, 2006;Ponette-Gonzalez et al., 2014). The rapid growth of various megacities located in the high Andes will further exacerbate the demand for water resources in the near future, from drinking water to A. Molina et al.: Multidecadal change in streamflow associated with anthropogenic disturbances water for sanitation, irrigation and agriculture, mining operations, and hydropower production. Changes in freshwater flow regimes are predicted to lead to future water scarcity (Bradley et al., 2006) as a result of the combined effect of climate change and variability (Bathurst et al., 2011;Urrutia and Vuille, 2009), and direct anthropogenic impact (Harden, 2006).
Both change and variability in local climatic conditions produce changes in hydrological conditions (Poveda and Mesa, 1997;Restrepo and Kjerfve, 2000;Poveda et al., 2011). In the tropical Andes, the temporal variability of precipitation is strongly related to oceanic and atmospheric conditions over the Pacific Ocean and Amazon Basin (Vuille et al., 2000;Marengo et al., 2004). The impact of El Niño-Southern Oscillation (ENSO) is clearly noticeable on the western escarpment of the Andes in Ecuador and northern Peru (Tapley and Waylen, 1990;Rossel, 1997;Vuille, 2013), and decreases with altitude as the steep, high-altitude topography of the Andean range creates distinct microclimates (Mora and Willems, 2012).
Direct anthropogenic impact associated with land cover change is rapidly transforming the hydrological functioning of tropical Andean ecosystems (Vanacker et al., 2003;Farley et al., 2004;Molina et al., 2012;Harden et al., 2013). The hydrological response is diverse, as changes in vegetation affect various components of the hydrological cycle including evapotranspiration (Nosetto et al., 2005), infiltration (Molina et al., 2007), and surface runoff (Bathurst et al., 2011;Restrepo et al., 2015). The clearance of native forest for arable and grazing land induces rapid changes in soil physical properties reducing soil infiltration capacity (Bosch and Hewlett, 1982;Molina et al., 2007) and increasing surface runoff as a result of soil compaction and reduced evapotranspiration (Ruprecht and Schofield, 1989). As a consequence, the conversion of native forests to agricultural land often results in an increase of the annual water yield but a reduction of the low flows (Bruijnzeel, 1990;Andréassian, 2004). In contrast, afforestation and/or reforestation of grasslands and arable lands lead to a reduction in soil moisture and total water yield as a result of greater canopy interception and evapotranspiration (Bruijnzeel, 2004;Scott et al., 2005;Farley et al., 2005;Buytaert et al., 2007).
In tropical Andean ecosystems characterized by large inter-and intra-annual variability in hydrometeorological conditions, little is known about the relative importance of climate change and direct anthropogenic perturbations on streamflow. At large spatial scales (> 100 km 2 ), the patterns of land cover change are notoriously dynamic, both in space and time, and are commonly associated with climatic and altitudinal gradients. In this paper, we assess multidecadal change in freshwater provision based on long time series of hydrometeorological data and land cover reconstructions. Given the strong temporal variability of precipitation and streamflow data related to ENSO, we use the Hilbert-Huang transformation to detrend the time series of streamflow and precipitation data. The adaptive data analysis is based on empirical mode decomposition techniques that are appropriate for nonlinear and nonstationary time series data (Huang et al., 1998). After empirical mode decomposition, the remaining long-term trends in streamflow and precipitation are contrasted to the observed patterns of land cover change.
The study is realized in an exceptional setting, the Pangor catchment (c. 282 km 2 ) in the Ecuadorian Andes. Situated on the western escarpment of the Ecuadorian Andes, the area is particularly affected by El Niño-Southern Oscillation cycles (Rossel, 1997). Land cover change resulted in a net loss of native forests and grasslands by about 20 % of the total catchment area between 1963 and 2009 (Balthazar et al., 2015). The main objective of this paper is to quantify the potential long-term effect of land cover change on streamflow in the tropical Andes. By analyzing multidecadal time series of hydrometeorological data, we specifically tested the relative sensitivity of streamflow to changes in land cover and precipitation.

Regional setting
The Pangor catchment (78 • 50 -79 • 01 W, 1 • 43 -1 • 58 S) is located 200 km southwest of the capital of Ecuador, Quito (Fig. 1). The catchment has pronounced relief with elevation ranging between 1434 and 4333 m a.s.l. over a distance of less than 30 km. Slopes are typically steep with slope gradients around 55 %, but steeper in the dissected river valleys. The climate can be described as equatorial mesothermic semi-humid to humid (Pourrut, 1994), with precipitation and temperature increasing strongly with altitude. Mean annual precipitation at the J. de Velasco station (3100 m a.s.l.; Fig. 1) is about 1400 mm , with high inter-annual variability ranging from 475 (2002) to 3700 mm (1994); whereas at the Chimbo DJ Pangor station (1450 m a.s.l.) annual precipitation is only about 1000 mm (INAMHI, 2009).
The underlying geology consists of volcanic and metasedimentary rocks of Cretaceous to early Tertiary age, with remnants of recent volcanic deposits at higher elevations. Soils have been classified as Andisols, Histosols, and Mollisols following the USDA (US Department of Agriculture) soil taxonomy (Gonzalez Artieda et al., 1986) and are characterized by a remarkably high water-holding capacity and soil organic matter content when undisturbed (Podwojewski et al., 2002). The landscape pattern now reflects several decades of land cover change. At mid and low altitudes, a complex patchwork of small agricultural plots, remnants of subalpine cloud forest, and patches of abandoned land with regeneration of natural shrub vegetation can be observed. Smallholder farming is the dominant agricultural activity, and crop rotation is a common practice where annual crops are alternated with pasture. Crop species vary with altitude, with maize (Zea mays) grown in association with the common bean (Phaseolus vulgaris) at altitudes below 2600 m a.s.l., and potato (Solanum spp.), faba bean (Vicia faba), and cereals (Triticum spp. and Hordeum vulgare) at higher altitudes. Large patches of montane cloud forest are the only remnant of forests on steep slopes in areas with very low accessibility. Above the natural treeline, the páramo grasslands are dominant, but plantation forests with exotic tree species (Pinus radiata and Pinus patula) now cover extensive areas (Balthazar et al., 2015).

Land cover change detection
Land cover change for the period 1963-2009 was reconstructed based on panchromatic aerial photographs (IGM, Quito, Ecuador) and high-resolution Landsat TM (15 October 1991) and ETM+ (3 November 2001, and 6 September 2009) images. A full coverage of aerial photographs at the scale of 1/60000 was obtained for November 1963 and 1977, and land cover mapping was realized following the orthorectification procedure described by Molina et al. (2012). Three Landsat scenes (1991,2001,2009, from the same season) with 1T level of pre-processing were acquired from the USGS (U.S. Geological Survey) archive, and images were atmospherically and topographically corrected with ATCOR3 (Atmospheric and Topographic Correction; Balthazar et al., 2012). To support the definition of land cover classes, a WorldView II image of 2010 with a horizontal resolution of 0.5 m (PAN) and 2 m (MS) was used (Digital Globe).
A multisource data integration method developed by  was applied to reduce imprecision and inconsistency that may result from the comparison of heterogeneous data sets (Balthazar et al., 2015). Four land cover types were defined: AL, agricultural land dominated by pastures and annual crops; F, montane cloud and subalpine forests (including primary and secondary forests); P, páramo grasslands dominated by tussock grasses and dwarf shrubs; and PP, exotic forest plantations dominated by Pinus radiata and Pinus patula.
The accuracy of the land cover change analysis is in function of the errors on the individual land cover maps. Land cover maps for 1963 and 1977 were extracted by manual on-screen digitalization on high-resolution copies of aerial photographs. As such, the accuracy of the land cover maps mainly depends on the horizontal positional accuracy of the orthorectified photographs and is systematically below the spatial resolution (30 m) of the aggregated land cover maps (Vanacker et al., 2003;Guns et al., 2014). A thorough validation of the land cover classification was realized for the 2001 satellite-derived map , based on a stratified sampling of 300 points for which the reference class was identified on very high-resolution aerial photographs. The error matrix reveals an overall classification accuracy of 94 % . Given the temporal dependence of land cover time series, the uncertainty on the amount of bi-temporal land cover change is not the limiting factor in our analysis.

Time series of precipitation and streamflow data
Hydrometeorological data were obtained from the National Institute of Meteorology and Hydrology of Ecuador (IN-AMHI). Time series of daily precipitation records are available for four meteorological stations located in or close to the Pangor catchment (J. De Velasco, Chimbo DJ Pangor, Pallatanga and Cañi-Limbe); daily streamflow data was obtained from the Pangor AJ Chimbo gauging station (Fig. 1). Hydrometeorological data of all stations is collected and recorded from manual readings, and data processing and quality control is realized by INAMHI. The time series are appropriate for studies of long-term trends in precipitation and streamflow, as they are of good quality (data gaps < 10 %) and cover a prolonged period of time (1974- Only during 1997 and 2003 is there a gap in the series of observed streamflow and precipitation data of more than 3 months. To obtain continuous time series data for EEMD (ensemble empirical mode decomposition) analysis, the multiple linear regression method described by Mora and Willems (2012) was used to fill data gaps. This method estimates correlation coefficients between all pairs of hydrometeorological stations, for each current and preceding month, and applies a multiple linear regression equation to predict missing flow or precipitation data. Given the low density of rain gauges in the Pangor catchment, we applied the regionalization method proposed by Mora and Willems (2012) to obtain catchment-wide or areal average precipitation depths. Based on data of altitude, vegetation pattern, and precipitation regime, four meteorological regions were delineated and the closest rain gauge station was assigned to each region ( Fig. 1). Additionally, an altitude correction factor was applied based on the observed positive linear relationship between mean annual precipitation and altitude (Mora and Willems, 2012). The areal average precipitation for the entire Pangor catchment was then calculated by summing the weighted precipitation (by surface area) of the four regions and dividing this value by the sum of the weights. The areal average daily precipitation depths, P d (mm), were aggregated into monthly data for the period 1974-2008.
The time series of streamflow data  is based on daily water stage readings at the Pangor AJ Chimbo gauging station (Fig. 1). Stage records (m) were converted into discharge records (m 3 s −1 ) using the stage-discharge rating curve developed by INAMHI. The daily discharge was then converted to daily water production of the catchment by multiplying the daily discharge (m 3 s −1 ) by 86 400 to convert the values to cubic meters. The equivalent water depth, WD d (mm), was then calculated by dividing the daily water production by the total catchment area (km 2 ) and multiplying by 1000 to convert the values to millimeters to allow direct comparison with precipitation records. The daily water depths were split into quick and slow flows using the hydrostatistical toolkit WETSPRO (Water Engineering Time Series PROcessing tool). This procedure is based on subflow separation techniques and applies a generalization of the original Chapman filter. The Chapman filter assumes exponential recession for the hydrological subflows and is derived from the general equation of a low-pass filter (Willems, 2009). To condense the time series data, the daily water depths were aggregated to monthly time steps.
To assess the consistency of the hydrometeorological data sets, two simple tests were carried out on the annual records of precipitation and equivalent water depth following Costa et al. (2003). Although it may be preferred to perform data quality assessment at the monthly timescale, the hydrometeorological data sets that were available for this study do not facilitate such task given the response time of the hydrological system. First, the mean annual evapotranspiration, ET yr , was estimated as the difference be- tween the mean areal average annual precipitation, P yr , and mean annual equivalent water depth, WD yr , for the time period 1974-2008, thereby assuming that long-term changes in soil water storage, S, can be neglected. This first estimate of ET yr of 1097 mm yr −1 or 3.0 mm day −1 (resulting from the subtraction of 559 mm streamflow yr −1 from 1656 mm precipitation yr −1 , 1974-2008) is consistent with published estimates of potential evaporation for this region of 1000-1100 mm yr −1 by the National Institute of Meteorology and Hydrology (INAMHI, 2009). Second, the annual evapotranspiration data, ET yr , are highly associated with the areal average annual precipitation depths, P yr (Fig. 2), which are characteristic for tropical Andean basins (Mora and Willems, 2012). These two assessments show that the quality of the hydrometeorological data is acceptable, so that they can reliably be used for analyses of inter-annual variability and long-term change. In this study, we refrain from discussing the seasonal or intra-annual variability in streamflow and precipitation.

Empirical mode decomposition (EMD)
Given the nature of hydrometeorological data in tropical Andean basins, which often display an abrupt pattern of amplitude and frequency modulation at different timescales, empirical mode decomposition, a technique developed by Huang et al. (1998), is an ideal method to extract physically meaningful signals in the nonlinear and nonstationary time series . The method separates nonlinear oscillatory patterns of higher frequencies from those of lower frequencies, until a constant or monotonic trend is ultimately obtained (Wu and Huang, 2004;Peel and McMahon, 2006). In contrast to more traditional time series analysis techniques, such as Fourier transformation and wavelet analysis, EMD is not based on linear and stationary assumptions (Huang and Wu, 2008). Huang et al. (1999) demonstrated that EMD is generally successful in retrieving the physically meaningful signals hidden in time series. In some Table 1. Proportions (in %) of the four land cover types for 1963, 1977, 1991, 2001, and 2009, and amount of land cover change for the period 1963-2009 as percentage of the catchment surface area (%). 1963 1977 1991 2001 2009 (1963-2009)  cases, however, the physically meaningful variability tends to spread into different intrinsic mode functions or IMFs (socalled mode-mixing problem as described by Huang and Wu, 2008). Wu et al. (2009) addressed this issue by performing noise-assisted data analysis which consists in building an ensemble of IMFs (ensemble EMD or EEMD) by introducing Gaussian white noise into the original signal before performing the EMD analysis techniques. The time series of monthly precipitation values and equivalent water depths covering the period 1974-2008 were decomposed into their intrinsic mode function (IMF i ) components and the residual or trend. Gaussian white noise (having an amplitude of 0.2 SD (standard deviation) of that of the original signal) was introduced in the original time series following Wu et al. (2009). The final IMFs and residual or trend were then computed as the average of their corresponding ensemble members, IMF i . Here, only the significant IMF i and their residuals or trends were averaged following the methods developed in Brisson et al. (2015). The significance is defined by the original trend having variability higher than the variability of a trend derived from a randomized time series. Here we propose a three-step method to derive the significance of the trend. Firstly, the monthly values of the observed time series, P and WD, were randomly distributed. The resulting time series features a variability similar to the observed one but without any meaningful trend. In total, 1000 random time series were generated from the empirical data. Secondly, following the process of the EEMD, a perturbation was added to the 1000 random time series. Finally, the trend in each random time series was derived using EMD. The original trend was defined to be significant if its variability is higher than the 99th percentile of the variability of the trends derived from the random perturbed time series.

Land cover type
Using EEMD analysis techniques, the time series of areal average monthly precipitation, streamflow, and baseflow were decomposed into six ensemble intrinsic mode function components plus their significant residual or trend. In total, 34 years of hydrometeorological data were processed.

Estimation of the long-term water balance
A budget approach was used to approximate the different components of the water cycle, including evaporation and transpiration (Bruijnzeel et al., 2006). First, the annual water balance for the entire catchment was approximated as P yr + HP yr = WD yr + ET yr + S, where P yr is the areal average precipitation (mm yr −1 ), HP yr is the horizontal rainfall and cloud interception (mm yr −1 ), WD yr is the equivalent water depth as derived from streamflow measurements (mm yr −1 ), ET yr is the evapotranspiration (mm yr −1 ), and S is the change in soil water storage in the catchment (mm yr −1 ). P yr and WD yr were established for the entire catchment, as no detailed hydrological information is available to further spatialize WD yr . Long-term changes in soil water storage, S, can be neglected, as soils are typically shallow on the western escarpment of the Andes so that deep infiltration is limited. Horizontal rainfall, HP yr , is also considered to be negligible for the catchment-wide water balance, as additional water input from the interception of cloud water and wind-driven rain is typically constrained to the belt of cloud forests (i.e., 11.6 % of the catchment area in 2009). We can then estimate the annual evapotranspiration as ET yr = P yr − WD yr .
Second, the long-term change in water balance was realized in the two ecosystems where major changes in land cover occurred (Table 1): the tropical montane cloud forest (defined as the landscape unit between 2200 and 3200 m a.s.l. originally covered by cloud forest), and páramo ecosystems (here defined as the entire landscape unit of high altitude above the continuous forest line, 3200 m a.s.l.). A partial water balance was computed for each of these ecosystems, following the methods described below. In the tropical montane cloud forest, the annual evapotranspiration in the cloud forests was estimated as ET yr = Et yr + Ei yr + Es yr , with Ei yr = P yr − T f yr + Sf yr (4) and Et yr = T f yr + Sf yr − WD yr , where Et yr is the transpiration loss or soil water uptake (mm yr −1 ), Ei yr the evaporation from wetted vegetation surface (interception evaporation) (mm yr −1 ), Es yr the evaporation from the soil surface (mm yr −1 , negligible under dense vegetation), T f yr throughfall and Sf yr stemflow (mm yr −1 ). Fleischbein et al. (2006) measured Tf and Sf for three catchments under the montane cloud forest in the southern Ecuadorian Andes. Results of Tf and Sf were on average 59 and 1 % of P , respectively. For all remaining land cover types, the annual evapotranspiration was estimated from the potential evaporation following the methods described in Allen et al. (1998): where K s is a water stress factor, K c is the crop coefficient, and ET o is the reference crop evapotranspiration. Here, we take the potential evaporation (PE) as determined by IN-AMHI (2009) as the reference crop evapotranspiration and do not apply an empirical correction factor to convert PE into ET o . The reason is twofold: (1) there is a good correspondence between the long-term ET (1097 mm) as estimated from the catchment water balance and PE (ranging between 1000 and 1100 mm, depending on altitude), and (2) the validation of evapo(transpi)ration data and models for the southern Ecuadorian Andes indicates that ET o determined by the Penmann-Monteith method corresponds within 10 % of the measured values of PE by INAMHI (Baculima et al., 1999). Crop coefficients for natural grass vegetation and agricultural land were established for the southern Ecuadorian Andes (Buytaert et al., 2006), and K c values were estimated at 0.42 for natural grassland and 0.95 for agricultural land. Pine plantations were attributed a K c value of 1 based on the crop coefficient established by Allen et al. (1998) for conifers. The water stress factor, K s , was set at 1, given that the annual PE (and ET o ) is lower than the long-term catchment average precipitation, P yr . Land cover data were used to estimate temporal changes in partial water balance over the period 1974-2008, as the main hydrological components are here parameterized based on land cover type. Land cover data for the Pangor catchment are discrete in time and provide information for 1963, 1977, 1991, 2001, and 2009 based on panchromatic aerial photographs and high-resolution Landsat images. A Markov chain model was used for the interpolation of the temporal land cover data for 1974 and 2008 using, respectively, land cover transition probabilities of 1963-1977-2009.

Land cover dynamics (1963-2009)
Three main land cover change trajectories can be distinguished: (1) expansion of agricultural land by an area equal to 14 % of the catchment area (or 39 km 2 ) in 46 years' time, (2) deforestation of native forests by an area equal to 11 % of the catchment area (or −31 km 2 ) corresponding to a mean rate of 67 ha yr −1 , and (3) afforestation with exotic species in recent years by about 5 % (or 15 km 2 ; Table 1, Fig. 3). Over the time period 1963-2009, about 50 % of the 64 km 2 of native forests was cleared and converted to agricultural land. Small forest remnants are now scattered over steep terrain and/or in poorly accessible sites at higher elevations. Deforestation rates were highest in the 1960s, 1970s, and 1980s with a net deforestation rate of 89 ha yr −1 , and slowed down to about 18 ha yr −1 for the period 1992-2009; similar to what was reported earlier for the Ecuadorian highlands (Vanacker et al., 2003;Guns and Vanacker, 2013). The pattern of afforestation stands in sharp contrast to the deforestation pattern (Fig. 3). Afforestation is mainly concentrated in the subalpine and alpine zones, and started in the early 1990s. About two-thirds of the total decrease in páramo grasslands (−23 km 2 ) results from exotic forest plantations, and only one-third from conversion to agricultural land.

Long-term trends in precipitation and streamflow (1974-2008)
The flow regime  largely mimics the yearly variation in precipitation, with maximum mean monthly streamflow in April (equivalent water depth of 86 mm) and low flow in September (25 mm). More than 60 % of the annual flow is concentrated in the period between February and June. Annual values of precipitation and streamflow reveal strong inter-annual variation (Fig. 4). Using EEMD, the times series of monthly precipitation values and equivalent water depths were decomposed into six ensemble intrinsic mode function (IMFs 1-6) components plus the significant residual or trend (Fig. 5). The EEMD detrending analysis shows that the precipitation and streamflow regime changed significantly over time (Fig. 6), as the residual trend is not flat. Based on the EEMD analysis , we can conclude that the observed changes in streamflow are not the result of long-term precipitation change, as the direction of the residual trends in streamflow and precipitation is opposite. Despite increased precipitation, there is a remarkable decrease in streamflow (Table 1). Over the period 1974-2008, the rate of change in streamflow and baseflow varied through time. Two periods of change can be distinguished based on the EEMD time series analysis with a transition occurring roughly at the beginning of the 1990s. Before the early 1990s, there is a notable decrease in monthly streamflow and baseflow while the precipitation amounts increase (Fig. 6). From the mid-1990s onwards, there is no systematic trend in streamflow, while the precipitation trend is still increasing (Fig. 6).
The mean daily water depth is about 0.4 mm lower in the period 1992-2008 compared to 1974-1991, despite an increase in the mean daily precipitation of 0.5 mm ( Table 2). The largest difference is observed for low water depths, with a decrease of the Q 95 and Q 90 by 77 and 75 %, respectively. The moderate water depths (Q 10 -Q 90 ) decreased by 24 %, and the highest ones (Q 1 ) only by 15 %. Results of a chi-square analysis indicate that changes in mean annual precipitation, streamflow, and evapotranspiration between the two periods are significant (p value < 0.005). Streamflow and evapotranspiration exhibit the largest change with −22 and +33 %, respectively. Interestingly, the magnitude of increase in estimated ET is three-fold greater than the increase in precipitation. When analyzing the monthly distribution of streamflow, it is clear that the largest decrease in streamflow is observed during the dry season (JJAS), followed by the first (JFMAM) and second (minor) rainy season (OND; Fig. 8a). Similarly, the estimated mean monthly baseflow is systematically lower (3-11 mm) during the most recent period (Fig. 8b). During the dry season (JJAS), about 60 % of the reduction in total flow can be attributed to the strong decrease in baseflow.

Changes in water balance for montane cloud forest and páramo ecosystems
Land cover dynamics observed in the Pangor catchment are characteristic for the tropical Andes, with conversion of montane cloud forests and páramo grasslands for agriculture and forestry. The time series of land cover data ) re- vealed shifts in land cover dynamics in the early 1990s. From 1963 to 1991, agricultural land increased rapidly at an annual rate of about 1 %, mainly as a result of deforestation of montane cloud forests at a rate of 2.08 % yr −1 , while few changes occurred in páramo grasslands. In the early 1990s, there is a shift from net deforestation to a net increase in forest cover, as a consequence of the deceleration of deforestation and strong increase in exotic-tree plantations (+15 km 2 ) in páramo grasslands (Table 1; Balthazar et al., 2015). Table 3 shows first-order estimates of changes in the partial water balance in montane cloud forest and páramo ecosystems over the period 1974-2008. The observed land cover changes in montane cloud forest and páramo ecosys- tems are estimated to have resulted in a net loss of annual water yield by 74 mm (or about 13 % of WD yr ) over the period 1974-2008. The development of 15 km 2 of pine plantation in páramo ecosystems is estimated to have increased transpiration losses by about 8.6 hm 3 or 31 mm (Table 3). Pine forests' water use is very high compared to native páramo vegetation as result of the large total leaf surface area and deep root systems (Buytaert et al., 2007), and it largely affects the soil wa- Table 3. Changes in catchment water budget over the period 1974-2008 in (a) montane cloud forest and (b) páramo ecosystems. ET past corresponds to the evaporative losses of the land units prior to land cover change, ET to the evaporative losses after land cover change, and ET to the overall change in evaporation due to land cover change during the period 1974-2008. The ET values are expressed in cubic hectometers to indicate the changes in partial water budgets for the land units undergoing land cover change. To allow direct comparison with the results from the long-term time series analyses at the catchment scale, ET is also expressed in millimeters by dividing the estimated water production of the land units (hm 3 ) by the total catchment area (km 2 ) and multiplied by 1000 to convert the values to millimeters. For cloud forest ecosystems, P corresponds to a hypothetical change in catchment average net precipitation, as a result of a reduction in the horizontal precipitation in cloud forests after clearance assuming that interception precipitation is 20 % of the measured precipitation.   ter storage and retention in organic-rich páramo soils (Farley et al., 2004). In addition, the conversion of ∼ 6 km 2 páramo grassland to agricultural lands is expected to have further increased the transpiration losses by 3.0 hm 3 or 11 mm (Table 3). Despite high solar radiation in the tropical Andes, the water use of native plants in páramo ecosystems is very low because of the evaporative characteristics of páramo grass species. The páramo grass tussock species can consist of up to 90 % of dead leaves, resulting in low ET values.
On the other hand, in montane cloud forests, the conversion of 50 % of the surface area of forest to agricultural land is expected to have had an impact on the annual ET through changes in plant evapotranspiration after forest removal ( Table 3). However, changes in cloud interception as a consequence of the removal of forests can largely outweigh the for- mer effects on the overall water yield, as suggested by Bruijnzeel et al. (2006) and Hamilton et al. (2008). With degradation or removal of cloud forests, the mass of moistureintercepting leaf surfaces, including epiphyte biomass on branches and stems, is lost, and horizontal precipitation from fog or cloud is consequently also reduced. In our calculations, we assume that the contribution of horizontal precipitation to the overall water budget is up to 20 % of ordinary rainfall based on Bruijnzeel (2004). As such, deforestation of native forests (by an area equal to 11 % of the catchment area) might have engendered a net loss of annual water yield, WD yr , by about 32 mm, mainly as a result of reduced atmospheric moisture input from fog or clouds (Table 3).

Soil hydrology following land cover conversions
Land cover conversions are often followed by a phase of intense soil degradation that further exacerbates the anthropogenic impact on surface hydrology . Soil erosion measurements based on fallout radionuclides for the Chimbo catchment (central Ecuadorian Andes) by Henry et al. (2013) clearly illustrate that soil erosion rates highly depend on land cover and management: erosion rates in páramo grasslands are estimated at 9 t ha −1 yr −1 , and are significantly higher in forest plantations, pastures, and croplands with erosion rates of, respectively, 21, 24, and 150 t ha −1 yr −1 . The latter values are similar to soil erosion estimates for highly degraded Andean environments in southern Ecuador (Molina et al., 2008;Vanacker et al., 2014). Accelerated soil erosion after land cover change has been shown to alter soil hydrological conditions, e.g., through a reduction of soil water infiltration rates and soil water retention capacity (Podwojewski et al., 2002;Molina et al., 2007). The analysis of flow duration curves established for the periods 1974-1991 and 1992-2008 indicate that the largest change in streamflow is observed for low flow depths (Fig. 7). About 60 % of the reduction in total flow results from decreasing baseflow (Fig. 8), which points to reduced soil water storage capacity in the Pangor Basin (Fig. 8). This observation suggests that land cover change not only affects the hydrological cycle directly through changes in transpiration or net precipitation but also indirectly through changes in soil hydrological conditions.

Conclusions
Land cover dynamics observed in the Pangor catchment are characteristic for the tropical Andes, with rapid deforestation of native forests and afforestation with exotic tree species in more recent decades. Given the nature of hydrometeorological data in tropical Andean basins, which often display an abrupt pattern of amplitude and frequency modulation at different timescales, EEMD is an ideal method to extract physically meaningful signals. The EEMD analysis shows that the observed changes in streamflow  are not the result of long-term change in measured precipitation. Despite increased precipitation, there is a remarkable decrease in streamflow that very likely results from anthropogenic disturbances that are associated with land cover change.
Over the period 1974-2008, the rate of change in streamflow varied through time. During a first phase (1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991), there is a notable decrease in monthly streamflow and baseflow while catchment average precipitation increases. The largest change in the catchment water balance is observed during this period of forest clearance at a rate of 2.08 % yr −1 . Model simulations using a partial water balance suggest that a 20 % reduction in atmospheric moisture input from fog or clouds as occult precipitation in cloud forests can contribute to a net loss of annual water yield by 7 % over the period 1974-1991. During a second phase (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), there is no systematic trend in monthly streamflow but a positive residual trend in catchment average precipitation, suggesting an increase in loss by evapotranspiration. At the same time, we observe a shift from net deforestation to net reforestation, as a consequence of the deceleration of deforestation and strong increase in pine plantations in páramo grasslands by an area equal to 5 % of the catchment area. The development of 15 km 2 of pine plantation is estimated to have increased transpiration losses by about 31 mm or 5 % of the annual water yield.
In conclusion, our analysis suggests that significant longterm change in streamflow can be associated with anthropogenic disturbances following land cover change. Land cover change not only affects the hydrological cycle directly through changes in transpiration or net precipitation but also indirectly through changes in soil hydrological conditions. Our observations point to the importance of land use planning, to minimize the potential impact of land cover change on freshwater flow regimes in the tropical Andes.