The impact of forest regeneration on streamflow in 12 mesoscale humid tropical catchments

Although regenerating forests make up an increasingly large portion of humid tropical landscapes, little is known of their water use and effects on streamflow (Q). Since the 1950s the island of Puerto Rico has experienced widespread abandonment of pastures and agricultural lands, followed by forest regeneration. This paper examines the possible impacts of these secondary forests on severalQ characteristics for 12 mesoscale catchments (23– 346 km2; mean precipitation 1720–3422 mm yr −1) with long (33–51 yr) and simultaneous records for Q, precipitation (P ), potential evaporation (PET), and land cover. A simple spatially-lumped, conceptual rainfall–runoff model that uses dailyP and PET time series as inputs (HBV-light) was used to simulateQ for each catchment. Annual time series of observed and simulated values of four Q characteristics were calculated. A least-squares trend was fitted through annual time series of the residual difference between observed and simulated time series of each Q characteristic. From this the total cumulative change ( Â) was calculated, representing the change in each Q characteristic after controlling for climate variability and water storage carry-over effects between years. Negative values of̂ A were found for most catchments andQ characteristics, suggesting enhanced actual evaporation overall following forest regeneration. However, correlations between changes in urban or forest area and values of Â were insignificant ( p ≥ 0.389) for allQ characteristics. This suggests there is no convincing evidence that changes in the chosenQ characteristics in these Puerto Rican catchments can be ascribed to changes in urban or forest area. The present results are in line with previous studies of mesoand macro-scale (sub-)tropical catchments, which generally found no significant change in Q that can be attributed to changes in forest cover. Possible explanations for the lack of a clear signal may include errors in the land cover, climate,Q, and/or catchment boundary data; changes in forest area occurring mainly in the less rainy lowlands; and heterogeneity in catchment response. Different results were obtained for different catchments, and using a smaller subset of catchments could have led to very different conclusions. This highlights the importance of including multiple catchments in land-cover impact analysis at the mesoscale.


Introduction
Tropical regions have experienced extensive changes in land use and land cover during the last few decades (Lepers et al., 2005). Continuously rising demands for crop land and timber have led to substantial deforestation in many regions (Drigo, 2005), and although the global tropical deforestation rate remains high at 13 Mha yr −1 (FAO, 2006), forest regrowth on abandoned agricultural land is increasing, particularly in Latin America (Aide and Grau, 2004;Hecht, 2010) and South East Asia (cf. Fox et al., 2000;Xu et al., H. E. Beck et al.: The impact of forest regeneration on streamflow 1999). Because these "secondary forests" account for approximately one-third of the total tropical forest area (Brown and Lugo, 1990;Hölscher et al., 2005), understanding the impact of forest regrowth on water yield is important for water resources management and planning purposes (Giambelluca, 2002;Bruijnzeel, 2004), and the development of viable "payments for ecosystem services" schemes (Landell-Mills and Porras, 2002;Lele, 2009). However, despite this recognized importance, little is known of the water use of secondary tropical forests, although there are indications of enhanced water use during the period of most active biomass accumulation (Giambelluca, 2002;Juhrbandt et al., 2004;Hölscher et al., 2005).
The relationship between forest cover and streamflow (Q) is subject to a long-standing and ongoing discussion (Andréassian, 2004), also in the tropics (Bruijnzeel, 2004;Calder, 2005). The influence of forest cover change on flooding is particularly contentious (e.g. FAO, 2005;Bradshaw et al., 2007;Van Dijk et al., 2009) whereas the effect of forestation on tropical dry-season flows is also under debate (Calder, 2005;Scott et al., 2005). The general contention is that the net effect of an increase in forest cover on dryseason flow depends on the "trade-off" between increases in Q due to enhanced soil water recharge on the one hand (as forestation generally increases soil macro-porosity and infiltration characteristics; Ilstedt et al., 2007;Bonell et al., 2010;Zimmermann et al., 2006Zimmermann et al., , 2010Hassler et al., 2011), and decreases in soil water reserves and Q on the other hand due to the higher water use of trees compared to crops, pasture, or scrubs (Bruijnzeel, 1989(Bruijnzeel, , 2004Jackson et al., 2005;Scott et al., 2005). Reviews of micro-scale (< 1 km 2 ) experimental catchment studies (e.g. Jackson et al., 2005;Brown et al., 2005) -mostly conducted outside the tropics and in non-degraded settings where soil infiltration characteristics are not likely to be improved substantially by forestation -suggest that the increase in vegetation water use is indeed more important, and thus an increase in forest cover commonly leads to a decrease in both total and dry-season Q. However, although direct experimental evidence of the "infiltration trade-off hypothesis" (Bruijnzeel, 1989;Bonell et al., 2010) is missing due to a lack of comprehensive studies, demonstrated reductions in amounts of headwater-or hillslope stormflow production after reforesting severely degraded land in various parts of the tropics (e.g. Chandler and Walter, 1998;Zhou et al., 2002;Zhang et al., 2004) should be large enough to overcome the associated increases in forest water use (Chandler, 2006;cf. Bruijnzeel, 2004;Scott et al., 2005;Zhou et al., 2010). Indeed, as long-term Q records for large, once degraded catchment areas are becoming available, evidence of improved baseflows (Q bf ) following large-scale land rehabilitation is beginning to be documented (Wilcox and Huang, 2010;Zhou et al., 2010).
The tropical island of Puerto Rico provides a unique opportunity to study the impacts of natural forest regeneration on Q at the mesocatchment scale (1-10 000 km 2 ), as high-quality and long-term hydro-climatic records and sequential land-cover data are available. Since the 1950s, Puerto Rico has seen widespread secondary forest regrowth on abandoned pastures, agricultural land (mostly sugar cane), and coffee plantations (Thomlinson et al., 1996;Aide et al., 2000;cf. Del Mar López et al., 1998). Although previous work on the relationship between land-cover change and Q using lumped meso-and macro-scale catchment data has experienced some difficulty demonstrating unequivocal results (e.g. Van Dijk et al., 2012), possibly stronger conclusions may be obtained by selecting catchments within similar regions (in this case central Puerto Rico; cf. Peña-Arancibia et al., 2010). On the whole, one would expect marked drops in total Q and Q bf during forest recovery in areas where the general extent of soil surface degradation before land abandonment was limited and soil structural characteristics (and thus infiltration opportunities) therefore remained relatively unaffected by forest regeneration (cf. Aide et al., 1996;Zou and Gonzalez, 1997). For catchments that experienced advanced soil degradation prior to agricultural abandonment, major declines in the volumes of both total Q and quickflow (Q qf ) would be expected during forest regrowth due to much improved infiltration and retention capacities (cf. Chandler and Walter, 1998;Zhou et al., 2002). The direction and magnitude of the change in Q bf will depend on the trade-off between the changes in vegetation water use and infiltration associated with forest regeneration (Bruijnzeel, 1989;Scott et al., 2005).
The following hypotheses are tested here: (1) there is a negative relationship between the area under regenerating forest and the change in total Q (i.e. Q qf + Q bf ); (2) Q qf shows a negative relationship with area under regenerating forest and a positive one with area under urbanization; and (3) depending on the trade-off between the changes in vegetation water use and infiltration associated with forest regrowth, Q bf shows either a negative, no, or a positive relationship with the area under regenerating forest. Specific objectives are to quantify the effects of forest regeneration and urbanization on total Q, Q qf , and Q qf .

Study area
Puerto Rico is an island with a tropical maritime climate, located in the north-eastern Caribbean occupying ∼ 8870 km 2 . The geology is dominated by the volcanic Cordillera Central, with a few major outcrops of plutonic rock (mostly granodiorite), and karstic limestones towards the far north and south (Olcott, 1999). Soils developed in the volcanic substrates are largely clayey Ultisols with a rapidly diminishing saturated hydraulic conductivity (K s ) with depth (Schellekens et al., 2004) whereas the granodiorites produce less clayey Ultisols with a less pronounced K s profile (Kurtz et al., 2011). Islandwide mean precipitation (P ) is ∼ 1700 mm yr −1 (Daly et al., 1994(Daly et al., , 2003. There is a moderate P seasonality, with the three driest months of the year (January-March) receiving on average ∼ 300 mm in total and the three wettest months of the year (September-November) receiving on average ∼ 775 mm in total (Daly et al., 1994(Daly et al., , 2003. The northern and eastern portions of the island receive ∼ 30 % more P due to the rising of the moisture-bearing trade winds against the slopes of the central mountain range (Calvesbert, 1970;García-Martinó et al., 1996;Daly et al., 2003).
During the second half of the 20th century, socioeconomic changes in Puerto Rico led to migration from (upland) rural areas and to (lowland) urbanization (Dietz, 1986). The associated abandonment of pastures and agricultural fields allowed secondary forests to develop over increasingly large areas as time progressed (Thomlinson et al., 1996;Grau et al., 2003;Helmer, 2004;Parés-Ramos et al., 2008;Fig. 1 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005  documented, both in terms of its expansion over time, and forest composition and structure Chinea and Helmer, 2003;Grau et al., 2003). Tree density typically reaches a peak between 25 and 35 yr after abandonment whereas species richness and forest structure resemble those of old-growth forest after ca. 40 yr of regeneration . Total actual evaporation (ET a ) of mature upland forest in the maritime tropical climate of Puerto Rico is high compared to tropical continental sites (Schellekens et al., 2000), mostly because of enhanced wetcanopy evaporation rates (Holwerda et al., 2006(Holwerda et al., , 2012, and the same may well apply to the island's secondary forests (cf. Giambelluca, 2002).
The locations of the 12 catchments examined here and the respective changes in land cover over the period 1951-2000 (see also Sect. 3.1) are shown in Fig. 1. The size of the catchments ranges from 23 to 346 km 2 (median size 42 km 2 ), mean annual P varies between 1720 and 3422 mm yr −1 (median value 2021 mm yr −1 ), and the length of simultaneous P , potential evaporation (PET), and Q records between 33 and 51 yr (median duration 44 yr; Table 1).

Land cover
Land-cover maps were obtained for 1951for , 1978for , 1991). The maps for 1951 (Brockmann, 1952;Kennaway and Helmer, 2007) and 1978 (Ramos and Lugo, 1994) were derived from aerial photography. Although the maps for 1951 and 1978 were rasterized at a resolution of ∼ 30 and ∼ 11 m, respectively, the actual mapping resolution used by the photo interpreters is estimated at ∼ 300 and ∼ 50 m, respectively. The 1991 and 2000 maps 1 (∼ 30 m resolution; Helmer and Ruefenacht, 2005) were derived from Landsat data using regression-tree modelling and histogram match-ing. All land-cover maps were reprojected to a common 0.0001 • (∼ 11 m) geographical grid by nearest-neighbour interpolation. To accommodate the different classification schemes used in the respective mapping exercises, each landcover class was assigned to a generalized class (Table 2). For each catchment the net changes in urban and forest areas from the start of the simultaneous P , PET, and Q records (Table 1) until 2000 were calculated by linear interpolation of urban-and forest-area time series.

Streamflow
All available daily Q records were downloaded from the US Geological Survey (USGS) National Water Information System 2 in December 2012, resulting in an initial dataset of 111 gauging sites. Catchment areas were derived for each site using the PCRaster software (Wesseling et al., 1996) and the USGS National Elevation Dataset (NED) digital elevation model (∼ 30-m resolution). The following criteria for inclusion here were applied: (1) the USGS published estimate of catchment area deviated by < 10 % from our computed catchment area; (2) the length of the Q record between 1950 and 2005 was > 30 yr; and (3) the catchment was not subject to flow regulation or affected by major anthropogenic water extraction. The latter was assessed using annual USGS Water-Data Reports 3 and a map of water supply intakes in the Luquillo Experimental Forest (Crook et al., 2007; an island-wide map of intakes is lacking). The final dataset comprised 12 catchments ( Fig. 1 and Table 1). For the conversion of measured discharge [feet 3 s −1 ] to areal mean Q [mm d −1 ] the computed catchment area was used. The following four Q characteristics were calculated, on an annual basis, to study changes in Q through time: (1) the annual 95th percentile (percent time not-exceeded) daily Q (Q p95 [mm d −1 ]; indicative of peak flows); (2) the annual mean Q (Q tot [mm yr −1 ]; indicative of total water yield); (3) the annual 5th percentile daily Q (Q p5 [mm d −1 ]; indicative of low   Calculated by dividing the diameter of a circle having the same area as the catchment, by the maximum length of the catchment (Schumm, 1956). c Record is defined as the period of simultaneous P , PET, and Q data. Several catchments have gaps in their record. d After Brockmann (1952) and Kennaway and Helmer (2007  For four catchments the daily Q strongly exceeded the daily P (see Sect. 4.1) on one or more days, indicating errors in the Q and/or P data. To prevent such errors from biasing the results, for some catchments parts of the Q record were excluded from the analysis. Specifically, for the Bauta catchment data for 1996-1998 were excluded, for the Grande de Loíza catchment data prior to 1961 were excluded, for the Fajardo catchment data for 1989 were excluded, and for the Inabón catchment data for 1975 were excluded.

Precipitation
Daily P data from the Global Historical Climatology Network-Daily (GHCN-D) database 4 (Gleason, 2002) and Table 3. HBV-light model parameter units, descriptions, and ranges used for the calibration. The K2 and MAXBAS parameters were not calibrated, but were set to 1 − k bf (cf . Table 1) and one, respectively. For the Valenciano, Canóvanas, Grande de Patillas, and Culebrinas catchments calibration ranges of 1.2-1.3, 0.7-0.8, 1.4-1.5, and 1.4-1.5, respectively, were used for PCORR. For the remaining eight catchments PCORR was set to one.  Daly et al., 1994Daly et al., , 2003 was used to ensure reliable long-term, elevation-corrected P for the catchments. The method used to obtain time series of daily P for each catchment is described in Sect. 4.1.

Minimum and maximum air temperature
Daily minimum and maximum air temperature (T min and T max , respectively [ • C]) were used to compute daily time series of PET for each catchment. T min and T max time series (> 20 yr) are available for 21 stations in Puerto Rico within the GHCN-D database (Gleason, 2002). Figure 2 presents the number of daily T min and T max observations between 1955 and 2010. Maps of mean annual T min and T max from the WorldClim dataset 7 (∼ 1 km resolution; Hijmans et al., 2005) were also used to ensure reliable long-term, elevationcorrected T min and T max for the catchments.  Figure 3 presents a flow diagram summarizing the various steps that were carried out to investigate whether changes in forest and/or urban area have influenced the observed Q characteristics. The various steps will be explained in detail hereafter.

Spatio-temporal interpolation and rescaling of climatic variables
Having reliable catchment-mean time series of the climatic variables (P , T min , and T max ) is important to prevent spurious trends in the simulated Q from influencing the results.
To calculate time series of the climatic variables from 1955 (marking the start of many records) to 2010 for the catchments an approach was developed that (1) removed the linear trend from the time series for each station and variable; (2) spatially interpolated the trend and daily-irregular components; and (3) reunited them on a per-pixel basis; after which (4) the time series were rescaled. More specifically, the following steps were carried out for each variable (where X denotes the variable in question): 1. T min and T max data were converted from • C to K. For each station with a record length > 20 yr, trends were calculated from annual mean X time series using the non-parametric Mann-Kendall statistical test (Mann, 1945;Kendall, 1975) with Sen's estimate of slope (Sen, 1968).
2. For each station daily X time series were de-trended and divided by the station mean such that the new mean is unity.
3. Daily maps with a resolution of 0.02 • (∼ 2 km) were computed from time series of X trend (having a constant value for each station) and from time series of de-trended unity X values using spatial interpolation (see next paragraph).
4. On a per-pixel basis the cumulative integral of the trend was calculated and offset such that the mean is unity, resulting in time series of net X change.
5. By multiplying the time series of net X change by time series of de-trended unity X values, time series of unity X values were obtained.
6. The time series of unity X values were multiplied by a map of elevation-corrected long-term mean X (PRISM for P , and WorldClim for T min and T max ; see Sects. 3.3 and 3.4, respectively).
7. Finally, daily X time series were calculated for each catchment by averaging over all the pixels comprising each catchment, and the T min and T max time series were re-converted from K to • C.
We employed the computationally efficient inversedistance weighting (IDW; Shepard, 1968;Dirks et al., 1998) technique for the spatial interpolation. The IDW technique requires a value for the distance-decay parameter, which controls how much weight is given to nearby stations relative to stations further away. Although it is recognized that this parameter varies per location (e.g. Lu and Wong, 2007), following recommendations by Garcia et al. (2008), and to reduce computational time, here a constant value of 3 was assumed.
The present approach has two advantages over the customary approach of using the nearest station with the longest record. Firstly, information from all nearby stations with records > 20 yr is incorporated. Secondly, our approach ensures unbiased mean annual values by rescaling against elevation-corrected long-term means. On the other hand, there are two caveats. First of all, portions of a time series may originate from different stations, thereby introducing spurious signals if, for instance, they have different seasonal patterns. Furthermore, for P , step six of the approach merely changes the intensity of the series, and does not correct for storm events unsampled by the isolated "point" network of meteorological stations.

Potential evaporation
The present study used the empirical Hargreaves equation (Hargreaves and Samani, 1985) according to guidelines set by the UN Food and Agriculture Organization (FAO; Allen et al., 1998) to assess long-term changes in the evaporative situation of the study catchments. The Hargreaves equation (Hargreaves and Samani, 1985) reads: where R a is the extraterrestrial radiation [mm d −1 ]. R a was computed as described by Allen et al. (1998). The approach followed to obtain catchment-wide daily mean time series of T min and T max is outlined in Sect. 4.1. The choice for the empirical Hargreaves method was motivated on the one hand by the lack of available long-term data for the climatic variables required for more physically based methods such as the Penman-Monteith equation (Monteith, 1965), and on the other hand by the better availability of T min and T max data. The Hargreaves method has been evaluated successfully against PET based on the Penman-Monteith equation (Trajkovic, 2007) and various other equations (Lu et al., 2005;Sperna Weiland et al., 2012).

HBV-light model
The HBV-light model (Seibert, 2005) was used to simulate Q. HBV-light is a spatially-lumped, conceptual rainfallrunoff model based on the HBV model (Bergström, 1976). HBV-light runs at a daily time step, has two groundwater stores and one unsaturated-zone store, and uses daily time series of P , PET, and T as inputs. T was not relevant in the present case since it is only used to drive the snow model subroutine. Rainfall interception is not estimated explicitly in HBV-light but is implicit in the BETA, FC, and UZL parameters.  Wilk et al., 2001). For in-depth discussion of the model, see Seibert (2005). Table 3 briefly describes the model parameters and lists the calibration ranges used.
We were unable to close the water balance for several catchments after exhausting all possible parameter combinations, probably due to errors in the PRISM P map, Q data and/or catchment boundary data, (non-quantified) water extractions, and/or inter-basin groundwater transfers. Therefore, an additional parameter (PCORR) was introduced that scaled P as required to match observed and predicted Q. Note that HBVlight does not consider groundwater flow within or between catchments.
Model parameters for each catchment were calibrated using Latin hypercube sampling (LHS; McKay et al., 1979). LHS is a more efficient alternative (Yu et al., 2001) to the commonly used Monte Carlo technique (Metropolis, 1987;Beven, 1993;Seibert, 1999). Using LHS the parameter space (Table 3) is split up in n equal intervals. Values for the parameters are generated by sampling each interval just once in a random manner. The model is run n times with random combinations of the parameter values from each interval for each parameter. Here n = 30 000 model runs were used to ensure convergence of the performance criterion. The first 10 yr of the record (Table 1) were used as spin-up period to initialize the groundwater stores after which the model was run for the entire period (i.e. the first 10 yr were run twice). The split-sample procedure was used to calibrate the parameters against data from the second half of the period and validate the parameters against data from the first half of the period (Refsgaard, 1997).
Parameters of hydrological models are commonly calibrated using a composite of objective functions (i.e. a summary statistic incorporating several measures of performance; e.g. Madsen, 2000). Here, three different objective functions that strike a balance between accurate representations of all portions of the hydrograph were used. The first objective function represents the Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970;NS [−]). The NS is defined as where Q s and Q o represent 3-day mean simulated and observed Q, respectively [mm d −1 ], t is the time step [−], and the summation is over t = 1, 2 For each run, the objective functions were combined to form a single aggregate measure using the combined-rank method (Booij and Krol, 2010). The main analysis was conducted using the 30 "best" Q simulations. Parameter uncertainty (e.g. Beven, 1993;Seibert, 1997) was quantified as described in Sect. 4.4. It should be noted that the K2 parameter in HBV-light was not calibrated, but was set to 1 − k bf (where k bf is the baseflow recession constant listed in Table 1). In addition, the MAXBAS parameter was not calibrated, but was set to one, since the travel time for the catchments under study is likely to be less than a day (confirmed by exploratory calibration efforts in all catchments). The PCORR parameter was calibrated only for catchments that had absolute VE values > 20 % for the calibration period when PCORR was set to one and parameter combinations were exhausted. For these catchments, the range of PCORR values used for the calibration was an initial estimate ±0.05, where the initial estimate was calculated based on the initially obtained VE values.

Evaluation of the impacts of land-cover change on streamflows
For each catchment and Q characteristic, annual time series of the deviation between observations and simulations were calculated as where D are the deviation time series [%], Q obs and Q sim are annual time series of the observed and simulated Q characteristics, respectively, x is the year, and d = 1, 2, . . . , 30 are the simulations. D integrates the effects of P and PET variability and basin carry-over storage on the Q characteristics, thereby isolating the impact of other factors, notably land-cover change.
For each catchment and Q characteristic, the trend (termed a J [%]) in the annual time series of D prior to 2000, when land-cover data are available, was calculated using the non-parametric jackknife resampling technique (Quenouille, 1956;Tukey, 1958).
where L [−] is the length of the record until 2000 (cf. Table 1).Â can be interpreted as the cumulative change in annual observed values of the Q characteristic after accounting for the effects of climate variability and carry-over water storage on the Q characteristic.
The standard error associated withâ J (termedσ J [%]) consists of two parts. The first (σ s [%]) was calculated from the mean dispersion ofâ within the time series of D, and is indicative of observational errors in P , PET, and Q data, nonlinearities in the land-cover impact on Q, and/or possible defects in the HBV-light model structure.σ s was calculated using the jackknife technique as follows: where std refers to the standard deviation. The two parts were combined to yield the standard error ofâ J (σ J [%]) as follows: The standard error associated withÂ (σÂ [%]) is given bŷ To examine whether changes in forest and/or urban area influenced the four observed Q characteristics, correlation coefficients were calculated between amounts of change in forest or urban area per catchment and corresponding values ofÂ. Conventionally a significance level of 0.05 is applied, but this level was adjusted for the number of inferences made (four in the present study) to 0.01 using the Bonferroni procedure (Bland and Altman, 1995).

Changes in climatic variables
Using the new P dataset , which comprises the records from 70 stations, trends for the island ranging from −0.25 to +0.41 % yr −1 (mean value +0.02 % yr −1 ) were found (Fig. 4a). In addition, a strong north-south disparity was observed, with positive P trends mainly identified to the south of the Cordillera Central, and negative trends north of it (Fig. 4a). Likewise, the new PET dataset , based on 21 stations, shows similar magnitudes of change, with per-pixel values ranging from −0.30 to +0.17 % yr −1 (mean value −0.03 % yr −1 ), but with a less clear spatial pattern compared to P (Fig. 4b). Table 4 lists the calibrated parameter values for the HBVlight model plus objective function scores for the calibration and validation periods for each catchment. Median Nash-Sutcliffe (NS) values of 0.64 and 0.63 (median of 12 values, where each value represents the mean NS value of the 30 best simulations for each catchment) were obtained for the calibration and validation periods, respectively (Table 4). The median absolute VE values were 3.1 and 7.2 % (median of 12 values, where each value represents the mean VE value of the 30 best simulations for each catchment) for the calibration and validation periods, respectively (Table 4). Catchments F, G, I, and L required optimization of the PCORR parameter (cf. Table 4). Figure 5 shows the land-cover time series for each catchment. In all catchments forest and urban areas increased markedly between 1951 and 2000 at the expense of pasture and agricultural areas. During the period of simultaneous P , PET, and Q recording (Table 1) the net change in forest and urban areas for the study catchments ranged from +2 to +55 % (mean value +26 %) and from +2 to +11 % (mean value +7 %), respectively. The changes in pasture and agricultural areas ranged from −19 to +26 % (mean value −7 %) and from −63 to +9 % (mean value −26 %), respectively.

Changes in streamflow characteristics
For each catchment and Q characteristic, Fig. 6 shows annual time series of the difference between observed and simulated Q as expressed by D (Eq. 4) and the corresponding fitted trend lines, whereas Table 5 lists values of the total changeÂ (Eq. 5) and the corresponding errorσÂ (Eq. 9). Catchments with low NS values generally gave higherσÂ values (cf. Tables 4 and 5). In general D tot shows the least temporal variability (Fig. 6), which is reflected by lowσÂ values forÂ tot ( Table 5). The low temporal variability of Q tot is likely due to the low degree of sampling error associated with the calculation of annual mean Q. Most catchments show progressive decreases over time in their time series of D (Fig. 6), as indicated by the mostly negativeÂ values (Table 5), suggesting decreases in observed Q characteristics that are unrelated to carry-over effects of water storage and climate variability. The most pronounced decreases in Q tot were found for catchments B, C, and J, whereas a moderate increase was found for catchment D. Pronounced decreases in both Q characteristics related to low flows (Q p5 and Q dry ) were found for catchments A, C, F, G, H, J, K, and L, whereas a (moderate) increase was found only for catchment I. Clear decreases in the Q characteristic related to peak flows (Q p95 ) were found for catchments B, C, I, and J, whereas strong increases were found for catchments D, H, and L. Figure 7 shows regressions between the amount of urban and forest area change vs. the cumulative change in annual time series of D prior to 2000 (expressed byÂ in Eq. 5) for each Q characteristic. In spite of strong increases in forest and urban area, and pronounced changes in D over time for the Q characteristics of several catchments, all correlations were insignificant (p ≥ 0.389). Nonetheless, a weak (i.e. non-significant) positive relationship can be observed between changes in forest cover and changes in annual total streamflow Q tot , when excluding catchments C and J (which appear to be outliers; Fig. 7c).

HBV-light model performance
The HBV-light performance was good for both the calibration and validation periods ( Table 4), suggesting that the simulated Q for the catchments can be used with confidence. Note that strong land-cover effects would have deteriorated performance statistics for the validation period. Several catchments required optimization of the PCORR parameter (cf. Table 4), possibly due to biases in the PRISM P map, uncertainties in the Q and/or catchment boundary data, water extractions, inter-basin groundwater transfers (potentially exacerbated by karst), and/or water recycling (e.g. Ellison et al., 2012), which combined may cancel out or amplify  Table 6. Previous studies on tropical, subtropical, and warm temperate meso-(≥ 1 km 2 ) and macro-scale (≥ 10 000 km 2 ) catchments investigating the impact of changes in forest area on Q, listed alphabetically by country name. Only findings based on time-series analysis were listed. Present study included for completeness. In the "Q change" column "No change" also includes insignificant changes in flows. Where possible the reported change in Q was corrected for P variability. Refers to changes in peak flows only. f Bruijnzeel (2004) argues that these increases in Q should be ascribed to urbanization rather than tree planting. g Represents the net forest area change of native forest cover (−9 and +2 %, respectively) and forest plantations (+6 and +4 %, respectively). h No relationship was found between the degree of forest-cover change and the change in  Table 2) land-cover classification maps for the years 1951, 1978, 1991, and 2000. one another. However, the influence of the P scaling on the results is probably limited because the simulated Q was only used to control for climate and storage carry-over effects.

Impacts of land-cover change on streamflows
We were unable to support the three hypotheses because no significant relationships (p < 0.01) were found between the change in forest cover or urban area, and the change in the investigated Q characteristics across the 12 catchments (Fig. 7). This result agrees with previous meso-to macroscale catchment studies in the tropics, subtropics, and warm temperate regions (Table 6), which mostly failed to demonstrate a clear relationship between Q and change in forest area. Our use of multiple catchments proved important since many individual catchments showed pronounced changes in flow characteristics that might well have been attributed to land-cover change otherwise.
Nevertheless, a weak positive relationship was observed between the change in forest cover and the change in Q tot (Fig. 7c), tentatively suggesting that regeneration of forests leads to increases in Q in Puerto Rico. If true, this would imply that in the investigated catchments the increase in Q due to enhanced rainfall infiltration during forest regrowth overrides the decrease in Q associated with the greater water use of the aggrading forests (cf. Bruijnzeel, 1989;Scott et al., 2005). The soils beneath young secondary forests in Puerto Rico show higher bulk density than beneath mature secondary forests (Weaver et al., 1987;Lugo and Scatena, 1995) and are less structured (Lugo and Helmer, 2004). Presumably, this reflects soil compaction by livestock and cropping prior to abandonment and subsequent secondary forest regeneration. However, the general level of soil degradation under pasture in Puerto Rico (cf. Aide et al., 1996) is probably not sufficient to cause widespread occurrence of overland flow, although surface erosion under sugar cane and coffee plantations has been reported to be rampant in parts of the island (e.g. Smith and Abruña, 1955;Del Mar López et al., 1998). Unfortunately, to date, no unequivocal case of a positive relationship between changes in Q and reforested degraded area has been reported in the tropics (cf. Table 6), although demonstrated decreases in the amounts of headwater-or hillslope stormflow generation after reforesting severely degraded land (Chandler and  1998; Zhou et al., 2002;Zhang et al., 2004;Sun et al., 2006) must be considered large enough to overcome the associated increases in forest water use (Chandler, 2006;cf. Bruijnzeel, 2004;Scott et al., 2005). Nevertheless, long-term Q data for large, once degraded but subsequently rehabilitated catchments in sub-humid Texas (Wilcox and Huang, 2010) and the humid Red Soils region of southern China (Zhou et al., 2010) have recently indicated gradually increased Q bf over prolonged periods of time following large-scale land rehabilitation and regreening, suggesting soil improvement to be the dominant factor in these cases. Further process-based work is required to substantiate this contention. Similarly, pending the results of hydrological process studies in Puerto Rico's regenerating forests (in particular, infiltration and soil water retention vs. forest water use) the presently obtained positive relationships between the change in forest cover and the change in Q tot may be spurious. The present investigation is confounded somewhat by modest urbanization occurring simultaneously with forest regeneration in some of the investigated catchments (Fig. 5). Since urbanization typically increases the area of impervi-ous surfaces, thereby enhancing the frequency and intensity of infiltration-excess overland flow, more Q qf is produced (e.g. Harto et al., 1998;DeWalle et al., 2000;Ziegler et al., 2004;Rijsdijk et al., 2007). This, if progressing beyond a critical threshold, can even result in reductions in Q bf (Van der Weert, 1994;Bruijnzeel, 1989Bruijnzeel, , 2004 on top of the reductions incurred already by the higher water use of the regenerating forest (Giambelluca, 2002). However, no significant relationships (p < 0.01) between the degree of urbanization and changes in any of the observed Q characteristics in our catchments were identified (Fig. 7). This is probably due to insufficient amounts of urbanization occurring (+2 to +11 %, mean value of +7 %; Fig. 5). Similarly, studies of the "flashiness" (sensu Baker et al., 2007) of the runoff behaviour of urbanized and forested catchments in Puerto Rico revealed no differences (Ramírez et al., 2009;Phillips and Scatena, 2010). However, Phillips and Scatena (2010) did detect significant differences (p < 0.05) in the frequency of stage change (cf. McMahon et al., 2003), possibly indicating locally enhanced Q qf due to urbanization.  A p95 (a-b), Q tot (c-d), Q p5 (e-f), and Q dry (g-h) values (Eq. 5; Table 5). Each letter in the scatter plots represents a different catchment. The change covers the start of the record (cf.  Wu et al. (2007) reported a change in Q tot of −25 % (corrected for changes in P using a single P station that was not included in the current study as it did not meet the necessary quality requirements), between 1973-1980 and 1988-1995 for the Río Fajardo catchment in north-eastern Puerto Rico (here identified as catchment H). Wu et al. (2007) attributed this negative change in Q to the small amount of forest regeneration occurring during that period (estimated at 8 % based on the interpolated time series shown in Fig. 5). For the same catchment, here a similar change in Q tot of −31 % (calculated using D tot ) was found between the respective periods. However, a much smaller estimated change of −5 % (standard deviation ± 4 %) was obtained between the respective periods when taking into account the complete record (48 yr), using the calculated trend in D tot prior to 2000 (â J ) and the corresponding standard error (σ J , Eq. 8). Thus, Wu et al. (2007) may have overestimated the reduction in Q tot for the Río Fajardo catchment because they only used part of the available record. Additionally, given that the Río Fajardo catchment appears to be an outlier in the present analysis (Fig. 7e, f) possibly due to non-quantified anthropogenic water extractions, any inferences for this catchment should be viewed with caution.

Field-based estimates of vegetation water use
To explore the possible changes in total water yield that may be associated with a conversion from shaded coffee, sugar cane or pasture (the dominant land uses in Puerto Rico in 1951, see Sect. 3.1) to (semi-)mature secondary forest (current situation), it is of interest to compare the typical amounts of total water use (ET a ) by the respective vegetation types. Field-based estimates of transpiration plus soil evaporation for shaded coffee and sugar cane in Puerto Rico range between 1.4-3.2 mm d −1 (Lin, 2010;Gutiérrez and Meinzer, 1994) and 2.5-4.6 mm d −1 (Vázquez, 1970;Fuhriman and Smith, 1951;Goyal and González-Fuentes, 1989;cf. Harmsen, 2003), respectively. The mean soil water uptake of various well-watered (lowland) pastures in northern Puerto Rico was estimated at 2.8 mm d −1 (Van der Molen, 2002), whereas for mature upland (tabonuco) forests, values converge around 3.0 (± 0.1) mm d −1 (Van der Molen, 2002;Wu et al., 2006). Thus, despite the deeper root systems of the forests (cf. Nepstad et al., 1994), mean soil water use for the agricultural crops under consideration, pasture, and forests is quite similar, possibly due to the relatively rainy climate prevailing in Puerto Rico all year round (Calvesbert, 1970). To this, rainfall interception losses (higher for forest) should be added. The best estimates of rainfall interception for mature tabonuco forest (∼ 21 % of mean annual P ; Holwerda et al., 2006) translate to ca. 1.2 mm d −1 for an annual P of 2000 mm yr −1 (i.e. the approximate average rainfall for the 12 study catchments), vs. 0.9 mm d −1 for shaded coffee (Siles et al., 2010) and 0.2 mm d −1 for sugar cane (Leopoldo et al., 1981). Extrapolating the above-mentioned average values to a year, gives mean estimated ET a totals for shaded coffee, sugar cane, and pasture of, respectively, ca. 1170, 1355, and 1020 mm yr −1 vs. ca. 1515 mm yr −1 for mature forest. Thus, the full conversion from pasture or crop land to mature forest could potentially change Q tot by about −160 mm yr −1 (in the case of sugar cane) to −495 mm yr −1 (in the case of pasture).
For the catchments examined here a mean modelled change in Q tot of −86 (± 124) mm yr −1 was found (calculated from values ofÂ tot and Q tot listed in Tables 1 and 5, respectively). If it is assumed that this change was solely due to the mean change in forest cover of +26 %, then the +100 % change in forest cover would have resulted in a change in Q tot of about −331 (± 478) mm yr −1 . Although the modelled change in Q tot lies in between the above estimates based on local field studies of ET a , the very high standard deviation of the estimated modelled change in Q tot suggests that this agreement may be mere chance. Therefore, the question re-emerges as to why analysis of the impacts of land-cover change on Q in mesoscale tropical catchments does not confirm the findings of micro-scale experimental studies, which clearly demonstrate a relationship between change in forest area and change in Q (Bosch and Hewlett, 1982;Bruijnzeel, 1990;Sahin and Hall, 1996;Brown et al., 2005;Jackson et al., 2005).

Potential explanations for the lack of relationships
Upscaling in hydrology is a challenging issue that has received considerable attention (e.g. Peterson, 2000;Rodriguez-Iturbe, 2000;McDonnell et al., 2007;Blöschl and Montanari, 2010). For many studies the lack of a clear relationship between deforestation and change in Q can be explained by the rapid regeneration of tropical forests, where ET a quickly reverts back to its original level and possibly exceeds it for decades (Giambelluca, 2002;Bruijnzeel, 2004;Juhrbandt et al., 2004). However, this explanation does not apply to the present study, since the catchments were exploited as pastures or agricultural fields for a sustained period of time prior to their abandonment and subsequent forest regeneration. Another possible reason for the inconclusive results obtained by many studies is the use of seasonal or annual mean Q characteristics only, whereas it is generally preferred to use characteristics related to the frequency and magnitude of Q (Alila et al., 2009(Alila et al., , 2010, particularly when evaluating the change in peak flows. Finally, the failure to correct for climatic variability by many studies may also have led to inconclusive or even erroneous outcomes. The first among the potential explanations for the lack of relationship between land-cover change and Q characteristics in Puerto Rico is covariance between land cover and climate in space, due to landscape ecology and history (cf. Van Dijk et al., 2012), or in time, due to land-cover changes altering the climate (Van der Molen, 2002;Pielke et al., 2007). Several mesoscale atmospheric modelling studies have produced conflicting results regarding the existence of such a relationship for Puerto Rico (Van der Molen et al., 2006;Comarazamy and González, 2011). The contrasting P trends reported here for different parts of the island (Fig. 4a) are also unable to provide evidence since similar magnitudes of forest regeneration occur across the island (Fig. 1). Also, because observed P and PET time series were input to the HBV-light model, most of the covariance between land-cover change, and P and PET should have been accounted for in the present analysis.
A second potential explanation is that the land-cover, P , PET, Q, and catchment boundary data used here contain more uncertainty than those used in small, commonly well instrumented, experimental studies. Using data for 278 Australian catchments Van Dijk et al. (2012) showed that the introduction of noise to long-term means of observed PET, P , and Q reduced the likelihood of detecting a land-cover signal. Here, as different methodologies were employed to derive the land-cover maps for the different years, and because the maps had different spatial resolutions and classification schemes (see Sect. 3.1), it is conceivable that the landcover time series contained uncertainties. For P a higher station density is probably desirable, particularly in tropical areas where the uncertainty in estimated catchment-wide P is exacerbated by the predominantly convective character of P and the limited spatial extent of individual storm cells (Nieuwolt, 1977;Hastenrath, 1991). Likewise, the PET data should also be viewed with caution, since net radiation, wind speed, and relative humidity, other dominant factors controlling PET magnitude and trends (see McVicar et al., 2012, and references therein), were not explicitly included in its calculation (cf. Eq. 1). The error in the Puerto Rico Q data is estimated to be 3-6 % (Sauer and Meyer, 1992). The cumulative error in annual P , PET, and Q was quantified bŷ σÂ (Eq. 9). Even thoughσÂ does not account for the error in the land-cover data or for drift errors in the P , PET, and/or Q data during the study period, values ofσÂ already constitute a substantial portion of the variance inÂ values (the estimated total change in the observed Q) for all Q characteristics (cf. Table 5 and Fig. 7). This suggests that errors in land cover, P , PET, and Q data may have restricted the detection of a relationship between land-cover change and Q characteristics in Puerto Rico. Kundzewicz and Robson (2004) and Chappell and Tych (2012) also suggested that a high degree of observational error may mask the identification of Q change.
A third potential explanation is that the amount of forest area change in the investigated catchments is less than that of most small catchment experimental studies. Based on small catchment experiments the critical threshold value for forest area change beyond which changes in Q can be detected is generally assumed to be ca. 20 % (e.g. Bosch and Hewlett, 1982). Increases in forest area for the investigated catchments were 2 to 55 % (mean 26 %), with increases exceeding 20 % in six of the twelve catchments. However, forest regeneration progressed from the headwaters to the lowlands (Grau et al., 2003) and by the time the records started (cf . Table 1) most of the headwaters would have been reforested already. Since most Q is produced in these rainier headwater areas (García-Martinó et al., 1996) the observed overall change in Q may be far less than expected on the basis of lumped representations of P and forest area.
A fourth potential explanation relates to anthropogenic water extractions from the investigated catchments. Irrigation on agricultural fields from local wells may have increased ET a , thereby masking the effect of increasing forest cover on Q. However, irrigation of sugar cane and pineapple in Puerto Rico occurs only on the drier southern coastal plains, and most irrigation water originates from lakes outside the study catchments (Molina-Rivera and Gómez-Gómez, 2008). Therefore, irrigation effects can probably be excluded as a possible explanation. However, water withdrawals associated with urbanization may have confounded the results, possibly causing reductions in Q bf that are unrelated to changes in forest cover. Although catchments affected by major water extraction were excluded from the analysis (see Sect. 3.2), it is possible that some catchments contain undocumented water intakes.
A fifth potential explanation, proposed by Zhou et al. (2010) and applicable to humid tropical settings, is that ET a under such conditions is constrained by PET (i.e. energy limited; Calder, 1998), and therefore the expected increases in ET a due to forest regrowth are not evident in the Q record. However, the interception evaporation component of ET a can be well in excess of PET, particularly on mountainous maritime islands (Schellekens et al., 1999;Roberts et al., 2005;Holwerda et al., 2006Holwerda et al., , 2012McJannet et al., 2007;Giambelluca et al., 2009). Moreover, an insignificant correlation was obtained between forest cover change and mean Q during the dry season (January-March; Q dry ), when PET is not expected to be a limiting factor (Fig. 7g). Hence, this is not seen as a convincing explanation.
Finally, a sixth potential explanation concerns heterogeneity between the study catchments in terms of vegetation cover and land-use history prior to land abandonment (not accounted for by the semi-quantitative classification of the present study), morphology, geology, and soils (McDonnell et al., 2007), influencing the response of Q to changes in forest cover. In addition, the use of lumped values for the climate variables at the catchment scale may not be valid due to spatial heterogeneity of P -PET-Q relationships (e.g. García-Martinó et al., 1996;Blöschl et al., 2007;Donohue et al., 2011).

Conclusions
Although there was considerable variability in the change in Q across the 12 examined catchments, the correlations between changes in urban or forest area and changes in Q were insignificant (p ≥ 0.389) for all four Q characteristics. Consequently, there is little evidence to support the hypothesis that there is a relationship between the degree of secondary forest regeneration and the change in Q for the investigated catchments. Nevertheless, most catchments exhibited decreases in Q tot , which may be attributable to enhanced vegetation water use associated with forest regeneration.
The most likely reasons for the lack of relationship between land-cover change and Q characteristics in Puerto Rico include (1) errors in the land cover, climate, Q, and/or catchment boundary data used in the analysis; (2) generation of Q is mostly in the rainier headwater areas that were already forested at the start of Q observations whereas subsequent changes in forest area mainly occurred in the drier lowlands; and (3) heterogeneity in the catchment response. Overall, it seems that the hydrological impacts of forest regeneration in Puerto Rico are currently unpredictable at the mesoscale, and further advances in our understanding are hampered by the quality, availability, and record length of the climatic and flow data. Our findings highlight the importance of catchment-scale analyses using multiple catchments but at the same time confirm the need for additional process studies at all stages of forest regeneration, notably of vegetation water use (both rainfall interception and transpiration) and changes in lowland hillslope hydrological response.