Flood risk reduction and flow buffering as ecosystem services – Part 2: Land use and rainfall intensity effects in Southeast Asia

Watersheds buffer the temporal pattern of river flow relative to the temporal pattern of rainfall. This “ecosystem service” is inherent to geology and climate, but buffering also responds to human use and misuse of the landscape. Buffering can be part of management feedback loops if salient, credible and legitimate indicators are used. The flow persistence parameter Fp in a parsimonious recursive model of river flow (Part 1, van Noordwijk et al., 2017) couples the transmission of extreme rainfall events (1−Fp), to the annual base-flow fraction of a watershed (Fp). Here we compare Fp estimates from four meso-scale watersheds in Indonesia (Cidanau, Way Besai and Bialo) and Thailand (Mae Chaem), with varying climate, geology and land cover history, at a decadal timescale. The likely response in each of these four to variation in rainfall properties (including the maximum hourly rainfall intensity) and land cover (comparing scenarios with either more or less forest and tree cover than the current situation) was explored through a basic daily water-balance model, GenRiver. This model was calibrated for each site on existing data, before being used for alternative land cover and rainfall parameter settings. In both data and model runs, the wet-season (3-monthly) Fp values were consistently lower than dry-season values for all four sites. Across the four catchments Fp values decreased with increasing annual rainfall, but specific aspects of watersheds, such as the riparian swamp (peat soils) in Cidanau reduced effects of land use change in the upper watershed. Increasing the mean rainfall intensity (at constant monthly totals for rainfall) around the values considered typical for each landscape was predicted to cause a decrease in Fp values by between 0.047 (Bialo) and 0.261 (Mae Chaem). Sensitivity of Fp to changes in land use change plus changes in rainfall intensity depends on other characteristics of the watersheds, and generalisations made on the basis of one or two case studies may not hold, even within the same climatic zone. A wet-season Fp value above 0.7 was achievable in forest–agroforestry mosaic case studies. Inter-annual variability in Fp is large relative to effects of land cover change. Multiple (5–10) years of paired-plot data would generally be needed to reject nochange null hypotheses on the effects of land use change (degradation and restoration). Fp trends over time serve as a holistic scale-dependent performance indicator of degrading/recovering watershed health and can be tested for acceptability and acceptance in a wider social-ecological context.


Introduction
Inherent properties (geology, geomorphology) interact with climate and human modification of vegetation, soils, drainage and riparian wetlands in effectuating the degree of buffering that watersheds provide (Andréassian, 2004;Bruijnzeel, 2004). Buffering of river flow relative to the spacetime dynamics of rainfall is an ecosystem service, reducing the exposure of people living on geomorphological floodplains to high-flow events, and increasing predictability and river flow in dry periods (Joshi et al., 2004;Leimona et al., 2015). In the absence of any vegetation and with a sealed surface, river flow will directly respond to the spatial distribution of rainfall, with only the travel time to any point of specific interest influencing the temporal pattern of river flow. Any persistence or predictability of river flow in such a situation will reflect temporal autocorrelation of rainfall, beyond statistical predictability in seasonal rainfall patterns. On the other side of the spectrum, river flow can be constant every day, beyond the theoretical condition of constant rainfall, in a watershed that provides perfect buffering, by passing all water through groundwater pools that have sufficient storage capacity at any time during the year. Both infiltration-limited (Hortonian) and saturation-induced use of more rapid flow pathways (inter and overland flows) will reduce the flow persistence and make it, at least in part, dependent on rainfall events. Separating the effects of land cover (land use), engineering and rainfall on the actual flow patterns of rivers remains a considerable challenge (Ma et al., 2014;Verbist et al., 2010). It requires data, models and concepts that can serve as an effective boundary object in communication with stakeholders (Leimona et al., 2015;van Noordwijk et al., 2012. There is a long tradition in using forest cover as such a boundary object, but there is only a small amount of evidence supporting this (Tan-Soo et al., 2016;van Dijk et al., 2009;van Noordwijk et al., 2015a).
In Part 1 , we introduced a flow persistence parameter (F p ) that links the two, asymmetrical aspects of flow dynamics: translating rainfall excess into river flow, and gradually releasing water stored in the landscape. The direct link between these two aspects can be seen from Eq. (4) in Part 1 : where Q t and Q t−1 represent river flow on subsequent days, P tx the precipitation on day t (or preceding precipitation released as snowmelt on day t) and E tx the preceding evapotranspiration since the previous precipitation event, creating storage space in the soils of the watershed. The first term on the right-hand side of the equation represents the gradual release of stored water, causing a slow decline of flow as the pools feeding this flow are gradually depleted. The second term reflects the part of fresh additions of water are partitioned over immediate river flow and the increase of stocks from which water can be gradually released. The derivation of the link depended on the long-term water balance, and thus assumed that all out-and inflows are accounted for in the watershed.
Commonly used rainfall-runoff models (including the curve number approach and SWAT models) only focus on the second term of the above equation (Ponce et al., 1996;Gassman et al., 2007), without link to the first. Various empirical methods for deriving "base flow" are in use, but details of the calculation procedure matter. Results in Part 1  for a number of contrasting meso-scale watersheds in Southeast Asia suggested that inter-annual variation in F p within a given watershed correlates with both the Richards-Baker (R-B) flashiness index (Baker et al., 2004) and the base-flow fraction of annual river flow. However, the slope of these relationships varied between watersheds. Here, in Part 2 we will further analyse the F p results for these watersheds that were selected to represent variation in rainfall and land cover, and test the internal consistency of results based on historical data: two located in the humid and one in the subhumid tropics of Indonesia, and one in the unimodal subhumid tropics of northern Thailand.
After exploring the patterns of variation in F p estimates derived from actual river flow records, we will quantify the sensitivity of the F p metric to variations in rainfall intensity and its response, on a longer timescale to land cover change. To do so, we will use a model that uses basic waterbalance concepts: rainfall interception, infiltration, water use by vegetation, overland flow, interflow and groundwater release, to a spatially structured watershed where travel time from sub-watersheds to any point of interest modifies the predicted river flow. In the specific model used, land cover effects on soil conditions, interception and seasonal water use have been included. After testing whether F p values derived from model outputs match those based on empirical data where these exist, we rely on the basic logic of the model to make inference on the relative importance of modifying rainfall and land cover inputs. With the resulting temporal variation in calculated F p values, we consider the time frame at which observed shifts in F p can be attributed to factors other than chance (i.e. null hypotheses of random effects can be rejected with accepted chance of Type I errors).

GenRiver model for effects of land cover on river flow
The GenRiver model (van Noordwijk et al., 2011) is based on a simple water-balance concept with a daily time step and a flexible spatial subdivision of a watershed that influences the routing of water and employs spatially explicit rainfall. At patch level, vegetation influences interception, retention for subsequent evaporation and delayed transfer to the soil surface, as well as the seasonal demand for water. Vegetation (land cover) also influences soil porosity and infiltration, modifying the inherent soil properties. Water in the root zone is modelled separately for each land cover within a subcatchment, the groundwater stock is modelled at subcatchment level. The spatial structure of a watershed and the routing of surface flows influences the time delays to any specified point of interest, which normally includes the outflow of the catchment. Land cover change scenarios are interpolated annually between time series (measured or modelled) data. The model may use measured rainfall data, or use a rainfall generator that involves Markov chain temporal autocorrelation (rain persistence). As our data sources are mostly restricted to daily rainfall measurements and the infiltration model compares instantaneous rainfall to infiltration capacity, a stochastic rainfall intensity was applied at subcatchment level, driven by the mean as parameter and a standard deviation for a normal distribution (truncated at 3 standard deviations from the mean) proportional to it via a coefficient of variation as parameter. For the Mae Chaem site in northern Thailand, data by Dairaku et al. (2004) suggested a mean of less than 3 mm h −1 . For the three sites in Indonesia we used 30 mm h −1 , based on Kusumastuti et al. (2016). Appendix B provides further detail on the GenRiver model. The model itself, a manual and application case studies, are freely available (http://www.worldAgroforestry.org/ output/genriver-generic-river-model-river-flow; van Noordwijk et al., 2011). Table 1 and Fig. 1 provide summary characteristics and the location of river flow data used in four meso-scale watersheds for testing the F p algorithm and application of the Gen-River model. Figure 1 includes a water tower category in the agro-ecological zones; this is defined on the basis of a ratio of precipitation and potential evapotranspiration of more than 0.65, and a product of that ratio and relative elevation exceeding 0.277. As major parameters for the GenRiver model were not independently measured for the respective watersheds, we tuned (calibrated) the model by modifying parameters within a predetermined plausible range, and used correspondence with measured hydrograph as test criterion (Kobold et al., 2008). We used the Nash-Sutcliffe efficiency (NSE) parameter (target above 0.5) and bias (less than 25 %) as test criteria and targets. Meeting these performance targets (Moriasi et al., 2007), we accepted the adjusted models as basis for describing current conditions and exploring model sensitivity. The main site-specific parameter values are listed in Table 2 and (generic) land-cover-specific default parameters in Table 3. Table 4 describes the six scenarios of land use change that were evaluated in terms of their hydrological impacts. Further description on the associated land cover distribution for each scenario in the four different watersheds is depicted in Appendix C.

Bootstrapping to estimate the minimum observation
The bootstrap method (Efron and Tibshirani, 1986) is a resampling methods that is commonly used to generate "surrogate population" for the purpose of approximating the sampling distribution of a statistic. In this study, the bootstrap approach was used to estimate the minimum number of observation (or yearly data) required for a pair-wise comparison test between two time series of streamflow or discharge data (representing two scenarios of land use distributions) to be distinguishable from a null hypothesis of no effect. The pair-wise comparison test used was Kolmogorov-Smirnov test that is commonly used to test the distribution of dis-  (Zhang et al., 2006). We built a simple macro in R (R Core Team, 2017) that entails the following steps: i. bootstrap or resample with replacement 1000 times from both time series discharge data with sample size n; ii. apply the Kolmogorov-Smirnov test to each of the 1000 generated pair-wise discharge data, and record the p value; iii. perform (i) and (ii) for different size of n, ranging from 5 to 50; iv, tabulate the p value from the different sample size n, and determine the value of n when the p value reached equal to or less than 0.025 (or equal to the significance level of 5 %); the associated n represents the minimum number of observations required.
Appendix D provides an example of the macro in R used for this analysis.

Empirical data of flow persistence as basis for model parameterisation
Inter-annual variability of F p estimates derived for the four catchments ( Fig. 2) was of the order of 0.1 units, while the intra-annual variability between dry and rainy seasons was 0.1-0.2. For all years and locations, rainy-season F p values, with mixed flow pathways, were consistently below dryseason values, dominated by groundwater flows. If we can expect F p,i and F p,0 (see Eq. 8 in Part 1, van Noordwijk et al., 2017) to be approximately 0.5 and 0, this difference between wet and dry periods implies a 40 % contribution of interflow in the wet season, a 20 % contribution of overland flow or any combination of the two effects.
Overall the estimates from modelled and observed data are related with 16 % deviating more than 0.1 and 3 % more than 0.15 (Fig. 3). As the Moriasi et al. (2007) performance criteria for the hydrographs were met by the calibrated models for each site, we tentatively accept the model to be a basis for a sensitivity study of F p to modifications to land cover and/or rainfall.

Comparing F p effect of rainfall intensity and land cover change
A direct comparison of model sensitivity to changes in mean rainfall intensity and land use change scenarios is provided in Fig. 4. Varying the mean rainfall intensity over a factor 7 shifted the F p value by only 0.047 and 0.059 in the case of Bialo and Cidanau, respectively, but by 0.128 in Way Besai and 0.261 in Mae Chaem (Fig. 4a). The impact of the land use change scenarios on F p was the smallest in Cidanau (0.026), intermediate in Way Besai (0.048) and relatively large in Bialo and Mae Chaem, at 0.080 and 0.084, respectively (Fig. 4b). The order of F p across the land use change scenarios was mostly consistent between the watersheds, but the contrast between the reforestation and natural forest scenario was the largest in Mae Chaem and the small-est in Way Besai. In Cidanau, Way Besai and Mae Chaem, variations in rainfall were 2.2 to 3.1 times more effective than land use change in shifting F p , whereas in Bialo its relative effect was only 58 %. Apparently, the sensitivity to changes in land use change plus changes in rainfall intensity depends on other characteristics of the watersheds, and generalisations made on the basis of one or two case studies may not hold, even within the same climatic zone.

Further analysis of F p effects for scenarios of land cover change
Among the four watersheds there is consistency in that the "forest" scenario has the highest, and the "degraded lands" the lowest, F p value (Fig. 5), but there are remarkable differences as well; in Cidanau the inter-annual variation in F p    . Frequency distribution of expected difference in F p in "paired-plot" comparisons where land cover is the only variable; left panels: all scenarios compared to "reforestation"; right panel: all scenarios compared to degradation; graphs are based on a kernel density estimation (smoothing) approach. is clearly larger than land cover effects, while in the Way Besai the spread in land use scenarios is larger than interannual variability. In Cidanau a peat swamp between most of the catchment and the measuring point buffers most of land-cover-related variation in flow, but not the inter-annual variability. Considering the frequency distributions of F p values over a 20-year period, we see one watershed (Way Besai) where the forest stands out from all others, and one (Bialo) where the degraded lands are separate from the others. Given the degree of overlap of the frequency distributions, it is clear that multiple years of empirical observations will be needed before a change can be affirmed. Figure 5 shows the frequency distributions of expected effect sizes on F p of a comparison of any land cover with either forest or degraded lands. Table 5 translates this information to the number of years that a paired plot (in the absence of measurement error) would have to be maintained to reject a null hypothesis of no effect, at 5 % probability (p). As the frequency distributions of F p differences of paired catchments do not match a normal distribution, a Kolmogorov-Smirnov test can be used to assess the probability that a nodifference null hypothesis can yield the difference found. By bootstrapping within the years where simulations supported by observed rainfall data exist, we found for the Way Besai catchment, for example, that 20 years of data would be needed to assert (at p = 0.05) that the reforestation scenario differs from agroforestation, and 16 years that it differs from actual and 11 years that it differs from degrade. In practice, that means that empirical evidence that survives statis-  tical tests will not emerge, even though effects on watershed health are real. At process level the increase in "overland flow" in response to soil compaction due to land cover change has a clear and statistically significant relationship with decreasing F p values in all catchments (Fig. 6), but both year-toyear variation within a catchment and differences between catchments influence the results as well, leading to considerable spread in the bi-plot. Contrary to expectations, the disappearance of "interflow" by soil compaction is not reflected in measurable change in F p value. The temporal difference between overland and interflow (1 or a few days) gets easily blurred in the river response that integrates over multiple streams with variation in delivery times; the difference between overland-or interflow and base flow is much more pronounced. Apparently, according to our model, the high macroporosity of forest soils that allows for interflow, and may be the "sponge" effect attributed to forest, delays delivery to rivers by 1 or a few days, with little effect on the flow volumes at locations downstream where flow of multiple days accumulates. The difference between overland-or interflow and base flow in time-to-river of rainfall peaks is much more pronounced.
Tree cover has two contradicting effects on base flow; it reduces the surplus of rainfall over evapotranspiration (annual water yield) by increased evapotranspiration (especially where evergreen trees or trees with a large canopy interception are involved), but it potentially increases soil macroporosity that supports infiltration and interflow, with relatively little effect on water holding capacity measured as "field capacity" (after runoff and interflow have removed excess water). Figure 7 shows that the total volume of base flow differs more between sites and their rainfall pattern than it varies with tree cover. Between years total evapotranspiration and base flow totals are positively correlated, but for a given rainfall there is a trade-off. Overall these results support the conclusion that generic effects of deforestation on decreased flow persistence, and of (agro)/(re)-forestation on increased Table 5. Number of years of observations required to estimate flow persistence to reject the null hypothesis of "no land use effect", at p value = 0.05 using Kolmogorov-Smirnov test. The probability of the test statistic in the first significant number is provided between brackets and where the number of observations exceeds the time series available, results are given in italics. flow persistence are small relative to inter-annual variability due to specific rainfall patterns, and that it will be hard for any empirical data process to pickup such effects, even if they are qualitatively aligned with valid process-based models.

Discussion
In the discussion of Part 1 , the credibility questions on replicability of the F p metric and its sensitivity to details of rainfall pattern vs. land cover as potential causes of variation were seen as requiring case studies in a range of contexts. Although the four case studies in Southeast Asia presented here cannot be claimed to represent the global variation in catchment behaviour (with absence of a snowpack and its dynamics as an obvious element of flow buffering not included), the diversity of responses among these four already point to challenges for any generic interpretation of the degree of flow persistence that can be achieved under natural forest cover, as well as its response to land cover change.
Where Fig. 8 in Part 1 (van Noordwijk et al., 2017) explored the relationship in inter-annual variation between flashiness index and F p in the actual data for the four watersheds, we can now repeat the analysis for the modelled results for each scenario. Figure 8 presents two examples with, again, evidence that the flashiness index and F p are related, but with considerable variation between the watersheds and a lower slope for the Cidanau watershed with its downstream flow buffering.
The empirical data summarised here for (sub)humid tropical sites in Indonesia and Thailand show that values of F p above 0.9 are scarce in the case studies provided, but values above 0.8 were found, or inferred by the model, for forested landscapes. Agroforestry landscapes generally presented F p values above 0.7, while open-field agriculture or degraded soils led to F p values of 0.5 or lower. Due to differences in local context, it may not be feasible to relate typical F p values to the overall condition of a watershed, but temporal change in F p can indicate degradation or restoration if a location-specific reference can be found. The difference between wet-and dry-season F p can be further explored in this context. The dry-season F p value primarily reflects the underlying geology, with potential modification by engineering and operating rules of reservoirs, the wet-season F p is generally lower due to partial shifts to overland and interflow pathways. Where further uncertainty is introduced by the use of modelled rather than measured river flow, the lack of fit of models similar to the ones we used here would mean that scenario results are indicative of directions of change rather than a precision tool for fine-tuning combinations of engineering and land cover change as part of integrated watershed management.
The differences in relative response of the watersheds to changes in mean rainfall intensity and land cover change suggest that generalisations derived from one or a few case studies are to be interpreted cautiously. If land cover change would influence details of the rainfall generation process (arrow 10 in Fig. 1 of Part 1, ; e.g. through release of ice-nucleating bacteria, (Morris et al., 2014;van Noordwijk et al., 2015b;Ellison et al., 2017) this can easily dominate over effects via interception, transpiration and soil changes.
Our results indicate an intra-annual variability of F p values between wet and dry seasons of around 0.2 in the case studies, while inter-annual variability in either annual or seasonal F p was generally in the 0.1 range. The difference between observed and simulated flow data as basis for F p calculations was mostly less than 0.1. With current methods, it seems that effects of land cover change on flow persistence that shift the F p value by about 0.1 are the limit of what can be asserted from empirical data (with shifts of that order in a single year a warning sign rather than a firmly established change). When derived from observed river flow data, F p is suitable for monitoring change (degradation, restoration) and can be a serious candidate for monitoring performance in outcome-based ecosystem service management contracts. Choice of the part of the year for which F p changes are used as indicator may have to depend on the seasonal patterns of rainfall.
In view of our results, the lack of robust evidence in the literature of effects of change in forest and tree cover on flood occurrence may not be a surprise; effects are subtle and most data sets contain considerable variability. Yet, such effects are consistent with current process and scaling knowledge of watersheds.
In summarising findings on the F p metric, we can compare it with existing ones across the seven questions raised in Fig. 1 of Part 1 . Comparator metrics can be derived from various data sources, including the amount (and/or quality) of forest cover upstream, the fraction of flows that is technically controlled, direct records of river flow (over a short or longer time period), records of rainfall and/or models that combine landscape properties, climate and land cover. Tentative scoring for these metrics (Table 6) suggests that the F p metric is an efficient tool for data-scarce environments, as it indicates aspects of hydrographs that so far required multi-annual records of river flow.

Conclusions
Overall, our analysis suggests that the level of flow buffering achieved depends on both land cover (including its spatial configuration and effects on soil properties) and space-time patterns of rainfall (including maximum rainfall intensity as determinant of overland flow). Generalisations on dominant influence of either, derived from one or a few case studies are to be interpreted cautiously. If land cover change would Figure 8. Relationship between F p value and R-B Flashiness index across years in four Southeast Asian watersheds under a "natural forest" and "degrading" scenario, simulated with the GenRiver model. Table 6. Comparison of metrics at various points in the causal network (Fig. 2 of Part 1,  that can support watershed management and prevention of flood damage on the list of seven issues (I-VII) introduced in Fig. 1 Part 1 * .
Integrated (5-7) terrain Terrain-based (7A and + climate + land use + 5 in Fig. 2 of Part 1) Based on river flow characteristics (4 in Fig. 2   influence details of the rainfall generation process this can easily dominate over effects via interception, transpiration and soil changes. Multi-year data will generally be needed to attribute observed changes in flow buffering to degradation/restoration of watersheds, rather than specific rainfall events. With current methods, it seems that effects of land cover change on flow persistence that shift the F p value by about 0.1 are the limit of what can be asserted from empirical data, with shifts of that order in a single year a warning sign rather than a firmly established change. When derived from observed river flow data, F p is suitable for monitoring change (degradation, restoration) and can be a serious candi-date for monitoring performance in outcome-based ecosystem service management contracts. Watershed health is here characterised through the flow pattern it generates, leaving the attribution to land cover, rainfall pattern and engineering of that pattern and of changes in pattern to further locationspecific analysis, in the same way a symptom of a high body temperature can indicate health, but not diagnose the specific illness causing it. The data sets analysed so far did not indicate that the flow persistence at high flows differed from that at lower flows within the same season, but in other circumstances this may not be the case and further care may be needed to use F p val-ues beyond the measurement period in which they were derived. While a major strength of the F p method over existing procedures for parameterising curve number estimates, for example, is that the latter depends on scarce observations during extreme events and F p can be estimated for any part of the flow record, the reliability of F p estimates will still increase with the length of the observation period.
Further tests on the performance of the F p metric and its standard incorporation into the output modules of river flow and watershed management models will broaden the basis for interpreting the value ranges that can be expected for well-functioning watersheds in various conditions of climate, topography, soils, vegetation and engineering interventions. Such a broader empirical base could test the possible use of F p as a performance metric for watershed rehabilitation efforts.  Rainfall 1989-2009, source: BWS 1998-2008, source: 1998-2002, source: 1976 2000-2009, source: 1954-2003, source: 1976-1998, source: data Sulawesi and PUSAIR KTI f ICHARM g PU and PUSAIR Reference 1 2 3 4 of detailed report Figure B1. Overview of the GenRiver model.
Appendix B: Genriver model for effects of land cover on river flow The Generic River flow (GenRiver) model (van Noordwijk et al., 2011) is a simple hydrological model that simulates river flow based on water-balance concept with a daily time step and a flexible spatial subdivision of a watershed that influences the routing of water. The core of the GenRiver model is a "patch" level representation of a daily water balance, driven by local rainfall and modified by the land cover and land cover change and soil properties. The model starts accounting of rainfall or precipitation (P ) and traces the subsequent flows and storage in the landscape that can lead to either evapotranspiration (E), river flow (Q) or change in storage ( S) (Fig. B1): The model may use measured rainfall data, or use a rainfall generator that involves Markov chain temporal autocorrelation (rain persistence). The model can represent spatially explicit rainfall, with stochastic rainfall intensity (parameters RainIntensMean, RainIntensCoefVar in Table 2) and partial spatial correlation of daily rainfall between subcatchments. Canopy interception leads to direct evaporation of an amount of water controlled by the thickness of water film on the leaf area that depends on the land cover, and a delay of water reaching the soil surface (parameter RainMaxIntDripDur in Table 2). The effect of evaporation of intercepted water on other components of evapotranspiration is controlled by the InterceptEffectontrans parameter that in practice may depend on the time of day rainfall occurs and local climatic conditions such as wind speed).
At patch level, vegetation influences interception, retention for subsequent evaporation and delayed transfer to the soil surface, as well as the seasonal demand for water. Vegetation (land cover) also influences soil porosity and infiltration, modifying the inherent soil properties. Groundwater pool dynamics are represented at subcatchment rather than patch level, integrating over the land cover fractions within a subcatchment. The output of the model is river flow, which is aggregated from three types of streamflow: surface flow on the day of the rainfall event, interflow on the next day and base flow gradually declining over a period of time. The multiple subcatchments that make up the catchment as a whole can differ in basic soil properties, land cover fractions that affect interception, soil structure (infiltration rate) and seasonal pattern of water use by the vegetation. The subcatchment will also typically differ in "routing time" or in the time it takes the streams and river to reach any specified observation point (with default focus on the outflow from the catchment). The model itself (currently implemented in Stella plus Excel), a manual and application case studies are freely available (http://www.worldAgroforestry. org/output/genriver-generic-river-model-river-flow; van Noordwijk et al., 2011. Figure C1. Land use distribution of the various land use scenarios explored for the four watersheds (see Table 4).