Reforestation-induced changes of landscape composition and configuration modulate freshwater supply and flooding risk of tropical watersheds

Impact of changes in land cover and land use on hydrological service of tropical watersheds is one of the focal research tropics in both hydrology and Land Cover Land Use Changes (LCLUC). Land fragmentation is an important feature of LCLUC, however, its impact on hydrological service of tropical watershed is unclear despite a few theoretical frameworks. In this paper, we described a simulation study of eight tropical watersheds in Puerto Rico using the Soil Water Assessment Tool. Annual average stream discharge was derived according to the simulations with the land cover maps in 1977, 1991, and 2000. Annual big stream discharge with risks of flooding and severe soil erosion was computed as the sum of daily discharge greater than 95th percentile. The impacts of changes in land cover and fragmentation represented by perimeter-to-area ratio of land patches on annual average and big discharges were further analyzed by means of the linear mixed-effects model. Most mountainous watersheds were characterized by reforestation in 1977–1991 but slight deforestation in 1991–2000. Forest perimeter-to-area ratio was significantly correlated with covers of forest (correlation coefficient of -0.97), pasture (0.94), and urban (0.95). Thus forest fragmentation was reduced by reforestation but increased by deforestation. The annual average and big discharges were significantly reduced by forest cover and forest perimeter-to-area ratio. The enhanced edge effect by forest fragmentation may have incurred more effective interception of the subsurface flow by forest root system, and promoted forest transpiration, thus reduced stream flows. Land cover change plays more important roles in regulating the big discharges than altering the annual average discharges. Due to the negative correlation between forest cover and fragmentation, the decreased forest fragmentation accompanied with reforestation offsets the impact of reforestation on lessening freshwater supply and flooding risk.

Introduction leave passages for overland or subsurface flow generated by other vegetation types from upper slope [31], a situation similar to the effects of fragmented Tigger Bush [32]. The net effect of changes in landscape configuration is largely unknown for regions with diverse reforestationdeforestation intensities.
This study intends to explore how changes in landscape composition and configuration impact the hydrological service of tropical watersheds in the Caribbean region. We attempt to test two hypotheses: (1) Freshwater service of tropical watersheds is affected by landscape composition and fragmentation, in additional to the main control of rainfall; and (2) Big stream flow with flooding and erosion risks, defined by daily discharge greater than 95 th percentile, is more affected by variation of landscape composition and fragmentation than the average stream discharge.
We used the Soil Water Assessment Tool (SWAT) [33][34][35] to simulate the discharges of eight tropical watersheds in Puerto Rico using three land cover/use maps from 1977 to 2000. The modeling results were checked with the USGS long-term gauged stream discharges. The average discharge and the big flow were further analyzed to derive relationship between these variables and the variables of landscape composition and configuration.

Methods Watersheds
We chose eight tropical watersheds in Puerto Rico based on the following criteria: (1) the watersheds should have USGS gauged discharge observations lasting for at least the past 20 years until 2010; (2) due to the lack of long-term water use data, watersheds with reservoirs are not considered except the Rio Grande de Manati watershed; and (3) the watersheds should have at least one meteorological station within or nearby with 80% or more complete daily rainfall and temperature records. The Rio Grande de Manati is an important watershed because the river originates from the highest peak in Puerto Rico and delivers large amount of freshwater to the metropolitan areas. We intended to do a water discharge simulation for this watershed despite two reservoirs, i.e., El Guineo and Matrullas, within it. The names, spatial locations, USGS discharge gauge codes, and brief description of selected watersheds were illustrated in Fig 1 and Table 1. These watersheds are mostly the mountainous upper reaches of large watersheds in Puerto Rico. Three watersheds, i.e., Rio Espiritu Santo (ESS), Rio Blanco (BLC), and Rio Fajardo (FJD), originate from the eastern Luquillo Mountains with great rainfall and high forest cover. Rio Grande de Loiza near Caguas (LOZ) and Rio Gurabo (GRB) are the upper parts of the Rio Grande Loiza watershed which spans from the mountains to the northern coastal wetlands. Rio Grande de Manati at Ciales (MNT) and Rio Cibuco (CBC) are located in the Karst area with unknown leaking to the underground streams. Rio Culebrinas (CLR) is located in the northwestern corner (Fig 1).

Data preparation
Hydrological simulation was carried out by means of the Soil Water Assessment Tools (SWAT) (http://swat.tamu.edu) which needs geomorphology, soil, and land use information as inputs. Daily rainfall and temperature from the meteorological stations within or nearby the watershed were utilized to drive the model. We preprocessed the DEM at 10m resolution (www.gis.pr.gov), the soil map and soil properties from the Soil Survey Geographic database (SSURGO) [36], and the land cover/use maps in 1977, 1991, and 2000 [37]. The land cover/use maps of 1991 and 2000 were digital-classified based on the Landsat images, but the map of 1977 was visually interpreted from the aerial photos [37] and has much lower resolution than the Landsat derived maps.
Daily climate records were mostly from the NOAA stations for all the watersheds except the Rio Espiritu Santo, for which the long-term rainfall records from the El Verde field station of the University of Puerto Rico were applied. Missing temperature data was interpolated based on the long-term averages, but missing rainfall was scaled from that of the nearest station. The 'borrowed' rainfall assumed daily synchronization among stations thus may incur error in the simulation.

The SWAT model and calibration
SWAT simulates the hydrological processes of a watershed together with the vegetation growth [33-35, 38, 39]. The hydrological components in the model include rainfall, soil water, shallow  groundwater, evapotranspiration, and surface/subsurface flows. The modeled stream discharge is usually validated against the gauged discharge at the outlet of watershed. Nutrient dynamics as well as sediment/chemical transportation in the water and the soil are also considered. SWAT is neither a lumped nor a pixel-based model. The basic simulation unit is the hydrological response unit (HRU) defined by the overlay of sub-basins, soil types, land cover/use, and slope categories. Surface and subsurface flows are routed to streams and rivers in the watershed. Parameters of soil and vegetation are connected to the built-in SWAT database and can be manually adjusted. The DEM, land cover/use, and soil data are processed by the ArcS-WAT module during the watershed delineation and the HRU definition. We set 200 hectare as the minimum sub-basin area to limit the details in the watershed delineation except BLC and ESS, for which 20 ha was used. We used three slope ranges with equal area as slope categories, and 15% as the minimum sub-basin area in the HRU definition.
SWAT was calibrated for each watershed manually for the discharge gauged during 1996-2005 and the land cover map of 2000. The curve number (cn2) was adjusted for appropriate runoff ratio. Coefficient of plant uptake (EPCO), soil evaporation (ESCO) from deep soil layer, and available soil water content (SOL_AWC) were adjusted for evapotranspiration. Saturated soil hydraulic conductivity (SOL_K), threshold water depth in the shallow aquifer required for return flow (GWQMN), coefficient controlling upward movement from shallow aquifer to overlaying unsaturated zone (QW_REVAP), and threshold depth of shallow aquifer allowing REVAP or percolation to deep aquifer (REVAPMN) were adjusted for appropriate lateral flow, evaporation from the shallow aquifer, and percolation to the deep aquifer (SWAT 2012 Input / Output document). All the parameters were adjusted within 20% of the nominal values.
Simulation using observed meteorological records from the 1970s to the end of 2013 For each watershed, we did simulation based on each of the three land cover/use maps with the same set of climate and soil data. The outputs included monthly and daily stream discharges at the outlets of the watershed and the sub-basins, as well as rainfall and evapotranspiration. The monthly discharge simulated with the land cover/use of 2000 was compared to the USGS gauged discharge. The quality of simulation was assessed by calculating Pearson's correlation coefficient between the gauged and the simulated discharges from 1996 to 2013, and root mean square of error. The 95 th percentile of daily discharge was calculated from the simulated daily discharge with the land cover/use in 1977. Daily discharges greater than the percentile were defined as big stream flows which tend to cause flooding and severe erosion, and then cumulated for each year.
During the post simulation analysis, we grouped the land cover/use types into four categories, i.e., urban, agriculture, grassland/pasture, and forest, and calculated the land cover fraction and the average perimeter-to-area ratio of each category. Variations in land cover fraction and perimeter-to-area ratio represent changes in composition and configuration, respectively. We divided the long-term average and big discharges for each watershed by watershed area to give flow per watershed area. The average annual and big discharges per watershed area were then regressed on rainfall, land cover fraction, and perimeter-to-area ratio using the linear mixed-effects model with watershed as a grouping variable. Linear mixed-effects model is intended to analyze structured data with a fixed effect component and a random effect component [40]. The fixed effect component gives the parameters for the entire population, and the random effect component is in general for individual experimental unit. In this study, we intended to derive the relationship between hydrological service and changes in land cover in the fixed effect component. We also wanted to fit a random intercept for each watershed to account for difference in factors other than changes in climate and land cover among watersheds. Because the resolution of the land cover/use map in 1977 is much coarser than that in 1991 and 2000, and calculation of fragmentation index is highly sensitive to the resolution [20], the application of linear mixed-effects model was carried out in two steps, i.e., the regressions of simulations using maps of 1977 and 1991 without, but those of 1991 and 2000 with, the perimeter-to-area ratios as explanatory variables. The free statistical software R [41] was used throughout the data preparation and the post simulation analyses. Functions of the LME in NLME package and the LMER in LME4 package were used. Backward stepwise selection of the independent variables was done manually to minimize the Akaike Information Criterion (AIC) for the fixed-effect component of the model.

Comparison of the simulated watershed discharges to the USGS observations
The simulated monthly average discharge rate in m 3 s -1 using the land cover/use map of 2000 was compared to the corresponding gauged discharge during 1996-2013 (Figs 2 and 3) to assess the simulation quality. Pearson's correlation coefficients between the observed and the simulated discharges ranged from 0.81 to 0.94, and root mean square of error varied from 0.5 to 5.3. The error increased with the amplitude of temporal variation (maximum discharge), and may be caused by the missing rainfall and temperature sequences and the unknown water use from the aquifer by local people. The unknown water use from the two mountain reservoirs may be responsible for the big error in the Rio Grande de Manati watershed. Land cover changes of the eight watersheds All watersheds except Rio Blanco had forest cover increased in 1977-1991 (Table 2). Except Rio Cibuco watershed which kept the reforestation, most watersheds experienced slight deforestation or stable forest cover from 1991 to 2000. The Rio Culebrinas watershed had the strongest reforestation with forest cover changing from 7% in 1977 to 41% in 1991. The Rio Blanco and the Rio Espiritu Santo watersheds maintained 90% forest cover due to the protection of the El Yunque National Forest. The average forest cover of all watersheds increased from 45% in 1977 to 58% in 1991 and then slightly decreased to 56% in 2000. Excluding the two protected watersheds (Rio Blanco and Rio Espiritu Santo), we found the average forest cover increased from 30% in 1977 to 46% in 1991, indicating more than 50% relative increase. Perimeter-to-area ratio was only computed for the maps of 1991 and 2000 (Table 2) due to the incompatible methodology used to derive the map in 1977. The forest perimeter-to-area ratio differed greatly among watersheds from 0.004 for Rio Espiritu Santo and Rio Blanco with 90% forest cover, to 0.03 m -1 for Rio Cibuco in 1991 and Rio Culebrinas in 2000. The perimeter-toarea ratio of forest in all watersheds maintained or increased from 1991 to 2000 except Rio Cibuco which was the only watershed keeping reforestation trend and thus had a decreased forest perimeter-to-area ratio.
Except the two watersheds in the El Yunque National Forest, the mountainous watersheds have relatively large pasture cover. The average grassland/pasture cover of the six watersheds outside the National Forest decreased from 44% in 1977 to 37% in 1991 and then increased to 39% in 2000. Decrease in pasture cover from 1977 to 1991 was associated with increase in forest cover for all the watersheds except Rio Culebrinas and Rio Cibuco for which pasture cover increased from the abandoned agriculture (Table 2). Pasture cover increased slightly from 1991 to 2000 for all the watersheds except Rio Cibuco and Rio Grande de Manati, and the increased pasture cover was associated with decreased forest cover. The average pasture perimeter-to-area ratio across all watersheds decreased slightly from 0.042 in 1991 to 0.04 in 2000. The increased pasture covers were all associated with reduced perimeter-to-area ratios, while decreased covers with enhanced perimeter-to-area ratios.
Most of the mountainous watersheds had agriculture cover less than 12% except Rio Culebrinas and Rio Grande de Manati, which had 67% and 24% agriculture in 1977 respectively. The agricultural cover of Rio Culebrinas decreased to about 3% and that of Rio Grande de Manati to 12%, giving rise to the forest, pasture, or urban development. However, there was a slight increase in agriculture during 1991-2000 for all the watersheds except Rio Gurabo and Rio Cibuco.
We found strong correlations among the land cover types (Fig 4). In particular, the Pearson's correlation between forest cover and pasture cover across all watersheds was -0.83, and that between forest cover and agricultural cover was -0.50, indicating the reforestation was mostly from pasture and abandoned agriculture.
Fragmentation was negatively correlated to the cover of the land cover type itself, but may be positively correlated with the cover of other land cover types (Fig 4). We found the forest perimeter-to-area ratio had a correlation coefficient of -0.97 with forest cover, and the pasture cover and its perimeter-to-area ratio had a correlation coefficient of -0.84. Therefore, reforestation reduced forest fragmentation. Forest perimeter-to-area ratio was positively correlated with covers of pasture and urban land, and the correlation coefficients were 0.94 and 0.95, respectively (Fig 4). Hence the fragmentation of forest may be caused by deforestation for Table 2. Landscape parameters and annual rainfall averaged over the 1981-2013 period for eight watersheds. A stands for land cover fraction of total watershed area, P is perimeter-to-area ratio, an index for shape irregularity and patchiness or fragmentation. Subscripts a, u, f, and g stand for agriculture, urban land, forest, and grassland/pasture, respectively.

WTSD
Year A a (%) pasture and urban expansion. Compared to forest, the perimeter-to-area ratio of pasture was positively correlated to forest cover with coefficient of 0.77, hence forest regrowth may increase the fragmentation of pasture. The perimeter-to-area ratio of pasture and that of forest also significantly correlated with a coefficient of -0.66. The abovementioned correlation coefficients were all significant at p<0.05 level.

Analyses of long-term averaged discharge and big steam flow as functions of land cover/use composition and configuration
For the simulations using the land cover/use maps of 1977 and 1991, we found the average annual discharge per watershed area F (Mt y -1 km -2 ) was significantly increased by the average annual rainfall, R (m), but reduced by the forest cover, C f (Eq 1).
Average discharge in the mountainous tropical watersheds was primarily controlled by rainfall variation which explained 91% of the total sum of regression squares. Changes in forest cover only explained 9%. The big stream flow per watershed area (F b ) was also increased by the rainfall but decreased by the cover of forest (Eq 2).
In contrast to the average discharge, the variation of forest cover explained 51% of the sum of regression squares, while the effects of rainfall largely decreased. When the two watersheds within the National Forest were excluded from the regressions, we found the coefficients were slightly altered (Eqs 3 and 4).
However, the forest cover explained 62% and 61% for the average and the big discharges, respectively.
For the simulations with the land cover/use maps of 1991 and 2000, the average stream discharge per watershed area was found as a function of rainfall and perimeter-to-area ratio of grassland/pasture (P g ) (Eq 5).
The average discharge was increased by rainfall, but reduced by the perimeter-to-area ratio of pasture which explained only 1% of the total sum of squares of the fixed effect model. The big flow per watershed area was found to be reduced by increased forest cover and its perimeter-to-area ratio (Eq 6), and the two variables accounted for 13% and 3% of the total sum of squares, respectively.
For the maps of 1991 and 2000, we also redid the regressions with exclusion of the two watersheds within the National Forest, and found different relationships (Eqs 7 and 8).
F 91 00 ¼ À 0:53 þ 0:87R À 0:31C f À 1:59P f ð7Þ The big discharge was not significantly affected by rainfall. Both the average and the big discharges were significantly decreased with forest cover and its perimeter-to-area ratio. Forest cover and its perimeter-to-area ratio explained 86% and 4% of the sum of squares for the average discharge, and 82% and 18% of that for the big discharge, respectively. All the coefficients in Eq (1) through Eq (8) had one-tailed p-values less than 0.05 according to the LME.

Discussion
The reforestation from abandoned agriculture and pasture in Puerto Rico is an important feature of the mountainous watersheds from 1977 to 1991. Increased forest cover in general reduces forest fragmentation. This conclusion is validated by the results at both individual watershed and multi-watersheds scales, and conforms to early findings [42]. In contrast to reforestation, deforestation due to urban sprawl and reclamation for agriculture increases forest fragmentation [11,22]. Human activities driven by the changes in policy and economy played the primary role in this transition [20,43].
Land cover changes largely impact stream discharges of the tropical mountainous watersheds. In 1970s-1990s, reforestation prevailed and forest cover explained around 60% of the total variation in the average and the big discharges per watershed area when the Rio Espiritu Santo and the Rio Blanco under national forest protection were excluded. When the two protected watersheds were included, the forest cover still explained 51% of the total sum of squares of the big discharges per watershed area. Our finding is comparable to previous studies [44] which concluded that landscape composition and configuration explained 64% of freshwater supply. Hence during the twenty years of fast reforestation, the increases in forest cover significantly reduced the stream discharge and the freshwater supply. On the other hand, the enhancement of forest cover also substantially reduced the big stream discharges, thus reduced the risk of flooding and landslides. The absence of the changes in the land cover type other than forest in the regression equation might be attributed to the prevailing signal of reforestation in this period. Land fragmentation indices were not applied in these regressions due to the incompatible resolution of the land cover map in 1977.
During late 90s to the new century, reforestation in Puerto Rico slowed down. The perimeter-to-area ratio of pasture appeared in the regression for average discharges despite its small contribution to the fixed component of the dependent variable (Eq 5). In the regression for big discharge and the regressions of average and big discharges without the two protected watersheds, forest cover accounted for up to 86%, and forest perimeter-to-area ratio for up to 18%, of the total variation of the dependent variables in the fixed models, in contrast to the previous conclusion that landscape configuration did not contribute significantly to the freshwater supply [44]. Decrease in forest cover during 1991-2000 significantly increased both average and big discharges for most watersheds. However, increase in forest fragmentation associated with the deforestation significantly reduced average and big discharges. The reason why forest cover did not appear in Eq 5 may have something to do with the great and stable forest cover of Rio Espiritu Santo and Rio Blanco, which may obscure signals of the changes in forest cover of other watersheds. Changes in forest cover and forest fragmentation have stronger impacts on big discharges than on average discharges, as shown by the greater percentages of explained variations in the regressions (Eqs 1-8). Forest has deeper root system and greater water budget than other land cover types, thus can buffer the extreme rainfall events and smoothen the stream flows. Our finding agrees with those found in the literature that reforestation reduces, but deforestation increases, the average discharge in the tropical watersheds [45,46]. Current literature showed equivocal consequence of changes in land configuration on hydrological service. Deforestation and forest fragmentation were found to increase hydrological connectivity [31] in the forest dominated landscape. On the other hand, forest fragmentation exposes more edges to solar radiation, so that water loss via transpiration is enhanced with fragmentation. More importantly, fragmented forest with more edges may have greater chances to intercept the overland and subsurface lateral flows along topographical slopes, so that forest can have more water to evaporate, than an aggregated forest. Hence the fragmented forest may also reduce the average stream discharges and the big flows. Our finding contradicts the argument that landscape fragmentation facilitates the delivery of ecosystem service [25,47].
Forest perimeter-to-area ratio is negatively correlated with forest cover. Reforestation is associated with reduced fragmentation, but deforestation is associated with increased forest fragmentation. When increased forest cover tends to reduce the stream discharges, the decreased forest fragmentation tends to enhance the stream discharges. On the other hand, when forest cover decreases during deforestation, the increased forest fragmentation tends to reduce the stream discharges. Thus changes in forestation fragmentation tend to offset the impacts of changes in forest cover in terms of regulation of freshwater supply and flooding risks.
The island of Puerto Rico has been under economic crisis in the past decade and the government decided to foster its agriculture sector. The recently issued land use plan calls for reserving at least 235,527 ha (582,000 acres) land with agriculture value which is mostly located in the central mountains (http://www.jp.pr.gov/). The land cover map of 2000 has already shown the deforestation for agriculture and pasture [20] despite the main trend of reforestation. The deforestation in the central mountains may continue to decrease the forest cover, but increase forest fragmentation. The net effect of these changes may increase the amount of freshwater supply, but on the other hand, raise the large discharges during episodic storm events.