Importance of landscape heterogeneity in sustaining hydrologic ecosystem services in an agricultural watershed

The sustainability of hydrologic ecosystem services (freshwater benefits to people generated by terrestrial ecosystems) is challenged by human modification of landscapes. However, the role of landscape heterogeneity in sustaining hydrologic services at scales relevant to landscape management decisions is poorly understood. In particular, the relative importance of landscape composition (type and proportion of land cover) and configuration (spatial arrangement of cover types) is unclear. We analyzed indicators of production of three hydrologic services (freshwater supply, surface and ground water quality) in 100 subwatersheds in an urbanizing agricultural landscape (Yahara Watershed, Wisconsin, USA) and asked: (1) How do landscape composition and configuration affect supply of hydrologic services (i.e., does spatial pattern matter)? (2) Are there opportunities for small changes in landscape pattern to produce large gains in hydrologic services? Landscape composition and configuration both affected supply of hydrologic services, but composition was consistently more important than configuration for all three services. Together landscape composition and configuration explained more variation in indicators of surface-water quality than in freshwater supply or groundwater quality (Nagelkerke/adjusted R: 86%, 64%, and 39%, respectively). Surface-water quality was negatively correlated with percent cropland and positively correlated with percent forest, grassland and wetland. In addition, surface-water quality was greater in subwatersheds with higher wetland patch density, disaggregated forest patches and lower contagion. Surface-water quality responded nonlinearly to percent cropland and wetland, with greater water quality where cropland covered below 60% and/or wetland above 6% of the subwatershed. Freshwater supply was negatively correlated with percent wetland and urban cover, and positively correlated with urban edge density. Groundwater quality was negatively correlated with percent cropland and grassland, and configuration variables were unimportant. Collectively, our study suggests that altering spatial arrangement of land cover will not be sufficient to enhance hydrologic services in an agricultural landscape. Rather, the relative abundance of land cover may need to change to improve hydrologic services. Targeting subwatersheds near the cropland or wetland thresholds may offer local opportunities to enhance surface-water quality with minimal land-cover change.


INTRODUCTION
Human domination of the biosphere is reshaping landscape heterogeneity that is critical for numerous ecological processes such as nutrient dynamics, energy flow and movement of organisms and materials (Turner et al. 2001, Foley et al. 2005, Kareiva et al. 2007).This, in turn, affects biodiversity and may undermine the capacity of landscapes to deliver ecosystem goods and services (benefits people obtain from nature) that are essential for human welfare (Millennium Ecosystem Assessment 2005, Tscharntke et al. 2005).Despite increasing understanding of the ecological basis of many ecosystem services, there remain significant knowledge gaps in the spatial ecology of ecosystem services (Nicholson et al. 2009, Mitchell et al. 2013).In particular, the links between landscape pattern and simultaneous provision of multiple ecosystem services, as well as the mechanisms underlying these relationships, are not well understood.This presents a significant challenge for management practices that attempt to alter landscape pattern to sustain and enhance ecosystem services in a changing world.
Hydrologic ecosystem services are defined as freshwater benefits to people generated by terrestrial ecosystems and include freshwater supply, water quality, flood mitigation, and water-related cultural services (Brauman et al. 2007).Hydrologic services are especially susceptible to landscape changes such as agricultural expansion or urbanization (Kepner et al. 2012).These changes can affect water quantity or quality by altering ecohydrological processes and introducing contaminants (Brauman et al. 2007, Kepner et al. 2012).For example, excess water extraction for intensive agricultural production or urban and suburban expansion can alter the water table and stream flows and lead to water scarcity (Foley et al. 2005, Carlisle et al. 2011).Excess nutrient additions to agricultural systems from commercial fertilizer and manure can trigger eutrophication and groundwater contamination (Carpenter et al. 1998).Expansion of urban impervious surface can increase ''flashiness'' of runoff from heavy rainfall events, resulting in more frequent and severe floods (Allan 2004).These landscape changes may interact with other drivers such as changing climate to challenge the sustainability of hydrologic services.
Understanding linkages between landscape pattern and the supply of ecosystem services is a looming research problem of import (Jones et al. 2013, Turner et al. 2013).A number of empirical studies have documented the role of landscape pattern for biodiversity-based ecosystem services such as pollination and pest regulation (Kennedy et al. 2013, Mitchell et al. 2014).However, most of these studies have focused on one service at a time, and few have examined multiple services, especially sets of hydrologic services.Different hydrologic services are likely to respond in varied ways to changes in landscape pattern, creating tradeoffs (i.e., negative relationships) or synergies (i.e., positive relationships) among ecosystem services.In consequence, an enhanced understanding of how landscape pattern affects the suite of hydrologic services and their interactions is needed for creating landscapes that can provide a balanced supply of hydrologic services.
Landscape pattern can affect hydrologic services directly by altering ecosystem processes such as hydrologic flow or lateral fluxes of nutrients, or indirectly by changing biotic communities.Previous research has contributed to understanding effects of land-use patterns on water quality and nutrient dynamics in different landscapes (Peterjohn and Correll 1984, Kesner and Meentemeyer 1989, Jones et al. 2001), and has elucidated the spatial scales over which these effects are manifest (Allan et al. 1997, Gergel et al. 1999, Strayer et al. 2003, Read et al. 2015).Nevertheless, these studies have placed greater emphasis on effects of landscape composition (type and proportion of land cover) than on effects of configuration (spatial arrangement of cover types).Empirical and theoretical evidence indicates that spatial configuration (e.g., connectivity, distribution, proximity or contagion of source and buffer ecosystems) may mediate the transport of water and nutrient across the landscape (Soranno et al. 1996, Weller et al. 1998, Gergel 2005, Alberti et al. 2007, Lee et al. 2009, Weller et al. 2011), thereby affecting hydrologic services.However, whether composition alone matters or if spatial configuration is also important for hydrologic services remains unresolved.Disentangling confounding effects between composition and configuration is challenging, and very few earlier studies of configuration effects controlled for effects of composition, producing conflicting or sometimes even misleading results (McGarigal andCushman 2002, Turner 2005).
Ecosystem functions and services are likely to respond nonlinearly to land-use changes (De-Fries et al. 2004, Turner 2005).An ecological threshold may occur at which there is an abrupt change in ecosystem quality, property or phenomenon, or where small changes in environmental drivers produce large ecosystem responses (Groffman et al. 2006).Prior studies have suggested that ecosystem thresholds can occur in a range of variables related to landscape pattern (Gergel et al. 2002, Turner 2005), and empirically estimated thresholds have been reported for landscape indicators such as impervious surface coverage (Paul and Meyer 2001), forest extent (Marzluff 2005) and mangrove area (Barbier et al. 2008) for stream health, bird species richness and coastal protection service, respectively.Understanding nonlinear responses to landscape pattern may suggest ecological leverage points where small management invest-ments can yield large benefits in ecosystem services (DeFries et al. 2004, Turner 2005, Bennett and Balvanera 2007) (Fig. 1).In particular, in areas where limited land-use/cover change is possible, identifying such nonlinearities can help target where and how to allocate these changes to maximize the production of ecosystem services (Gret-Regamey et al. 2014).However, detection of such nonlinearities remains challenging, and empirical evidence is still scant.
We studied the production and spatial distributions of indicators of three hydrologic ecosystem services (freshwater supply, surface-water and groundwater quality; Table 1) in the Yahara Watershed, Wisconsin (Fig. 2), an urbanizing agricultural landscape in the Midwest that typifies social-ecological stresses (e.g., climate change, urbanization, agricultural intensification, human population growth) on freshwater ecosystems and the services they provide (Carpenter et al. 2007).We analyzed effects of landscape composition and configuration on hydrologic services at the subwatershed scale, which is relevant for many land management decisions.Specifically, we asked: (1) How do landscape composition and configuration affect supply of Table 1.Biophysical indicators, median and range for three hydrologic ecosystem services quantified and mapped at the subwatershed scale for the Yahara Watershed, Wisconsin, USA.
hydrologic ecosystem services (i.e., does spatial pattern matter)?(2) Are there opportunities for small changes in landscape pattern to produce large gains in the provision of hydrologic services?We expected an overall strong influence of landscape composition, with groundwater recharge (the indicator for freshwater supply service) negatively correlated with percent urban cover, and water-quality services negatively correlated with the amount of agricultural and urban covers and positively correlated with natural land covers.We also expected landscape configuration to mediate supply of services dependent on lateral transfers of water and nutrients.Specifically, we hypothesized a strong edge effect on reducing surface flow rate, thus enhancing infiltration and groundwater recharge, whereas well-connected urban covers would increase runoff and reduce recharge.We also expected a strong positive edge effect of natural covers (e.g., forest or wetland) on facilitating nutrient absorption and retention.In addition, we expected isolated and less connected parcels of agricultural land (e.g., cropland or grassland) to reduce nutrient loss along flow pathways, as compared to the same number of parcels that were well connected.We also expected disaggregated and widely dispersed natural land covers to increase opportunities to trap nutrients and attenuate nutrient loss.At the v www.esajournals.orgsubwatershed level, we expected that subwatersheds with a finer-scale interspersion of landcover types (i.e., low contagion) would reduce nutrient loss by increasing interaction among different land covers.Please refer to Table 2 for the full list of hypotheses regarding relationships between landscape pattern variables and the indicators of each hydrologic service.

STUDY AREA
The Yahara Watershed in southern Wisconsin, USA (at 4686 0 N, 89824 0 W), is home to 372,000 people, drains 1,345 km 2 and contains five major lakes (Mendota, Monona, Wingra, Waubesa and Kegonsa; Fig. 2).The climate is continental, characterized by warm humid summers and cold winters with strong seasonal and inter-annual variability.Mean temperature ranges from 7.38C in January to 21.88C in July, with mean annual precipitation of 87.6 cm (1981-2010 climate normal; data available at www.ncdc.noaa.gov).Influenced by the last glaciation (;14,000 yr ago), the terrain is generally flat and marked by moraines, drumlins and shallow depressions (Auclair 1976).Soils are primarily composed of fertile Mollisols and Alfisols, which support high agricultural productivity (Glocker and Patzer 1978).
The Yahara Watershed is currently human dominated; most land is agricultural, but the region also includes densely populated urban and suburban areas, typical of many Midwestern landscapes (Fig. 2).The watershed has under-gone two major waves of landscape change since 1800 (Carpenter et al. 2007).Prior to European settlement (;1830s), this landscape was a mosaic of prairie, oak savanna, forest and wetland.Native vegetation was cleared for agriculture during the mid-1800s (Curtis 1959), and by 1870 most arable land was in production.The second wave became prominent in the mid-1900s as urban and suburban development expanded at the expense of formerly agricultural lands.At present, the watershed is about 60% agricultural (principally corn, soybeans and dairy farms), 25% urban (including Madison, Wisconsin, and adjoining towns), and includes scattered remnants of natural and semi-natural vegetation.The mix of different land covers makes for a heterogeneous landscape that interacts strongly with water.
The Yahara Watershed generates multiple hydrologic services that are important at local and regional scales.Freshwater is central to this region's cultural identity and wellbeing, but managing freshwater in this watershed has been remarkably challenging in the face of problems such as eutrophication, flooding, groundwater extraction and contamination.Eutrophication of the Yahara lakes has been an issue since the late-1800s, initially as a consequence of sewage effluent and erosion.In the late 1940s, eutrophication worsened due to excess nutrient inputsprimarily phosphorus and nitrogen-from increased sewage discharge, fertilizer and manure application, and agricultural and urban runoff (Carpenter et al. 1998, Lathrop 2007).Nutrient Note: Phosphorus loading is the inverse indicator of surface-water quality service.
v www.esajournals.orgenrichments also triggered groundwater pollution with excess nitrate leached through the soils.
Groundwater extraction, draining of wetlands, and increased runoff associated with expansion of impervious surfaces have altered hydrology, affected availability of usable clean waters, and increased lake-level variability (Soranno et al. 1996).Institutions for managing freshwaters have emerged since 1950s.Since 1970s, substantial efforts have been expended to curb freshwater problems at state and county levels (Wardropper et al. 2015), including wastewater diversion, biomanipulation, erosion control, storm water management, nutrient management, and wetland restoration (Carpenter et al. 2007, Lathrop 2007).Private groups and conservation organizations also have raised public awareness of freshwater issues and advocated for sustainable freshwater management.While many collective policy and practice efforts have been made to improve freshwater conditions, persistent changes in key drivers (e.g., land-use/cover change) and their long-term legacies still pose substantial challenges for the sustainability of hydrologic services.

METHODS
Our analysis unit was the subwatershed-a spatial scale at which landscape pattern is likely to affect ecohydrological processes and hydrologic services (Strayer et al. 2003).Second-order subwatersheds (n ¼ 100) were delineated based on 1:24,000 scale stream network, light detection and ranging (LiDAR) elevation, and a fieldchecked basin map from Dane County, Wisconsin.Subwatersheds were delineated using ArcHydro module in ArcGIS 10.0 and averaged 1,270 ha in area (standard deviation ¼ 565 ha; Appendix: Fig. B1).We averaged indicators of hydrologic services estimated at 30-m resolution (see Quantifying hydrologic ecosystem services) for each subwatershed as response variables, and calculated selected landscape metrics using 30-m land use/cover data (see Landscape pattern analysis) for each subwatershed as predictor variables, which were then analyzed (see Statistical analyses) to examine relationships between landscape pattern and hydrologic services.

Quantifying hydrologic ecosystem services
We used biophysical indicators to quantify three hydrologic services (Table 1) that are vital for the Yahara Watershed and for which spatial data and validated models were available.The selected indicators capture key ecohydrologic processes that underlie production or condition of each service and were mapped by Qiu and Turner (2013) at 30-m spatial resolution using empirical estimates and spatially explicit models for 2006-the most recent year for which data exist for all services (Appendix: Fig. C1).For the current study, we averaged each indicator by subwatershed (Fig. 3).Please refer to Qiu and Turner (2013) for full details on data sources, methods and accuracy assessment for each service.Below we provided a brief description of the approach for quantifying each service.
Freshwater supply (annual groundwater recharge, cm yr À1 ).-Groundwater recharge was used as the indicator for freshwater supply because this region depends on groundwater for industrial and municipal uses (e.g., domestic, irrigation, livestock and industrial relied on 100%, 97.5%, 90% and 84%, respectively on groundwater; Buchwald 2011).Hence, sufficient recharge to replenish aquifers is critical to sustain current and future human demands for freshwater in this watershed.Groundwater recharge was estimated using modified Thornthwaite-Mather Soil-Water-Balance (SWB) model (Dripps and Bradbury 2007) that was validated for this region.SWB is a physical-based quasi 3-dimensional model that accounts for processes of precipitation, evapotranspiration, interception, surface runoff, soil moisture storage and snowmelt.Groundwater recharge was calculated on a grid cell basis at a daily time step, with the following mass balance equation (Dripps and Bradbury 2007) We ran the model for a 3-yr period (2004)(2005)(2006) at 30-m resolution (Appendix: Fig. C1), with the first two years as ''spin-up'' of antecedent conditions that influence groundwater recharge for the focal year of 2006.We compared our recharge results with a recent 10-yr average estimates (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) for Dane County, and our estimates matched reasonably well with averaged groundwater recharge at the sub-basin levels v www.esajournals.org(Qiu and Turner 2013; see Appendix).
Groundwater quality (probability of groundwater nitrate concentration , 3.0 mg L À1 , unitless 0-1).-Nitrate is the most ubiquitous groundwater contaminant that has detrimental impacts on human health, and nitrate concentration is a useful indicator of groundwater quality (Spalding and Exner 1993).We used the probability of groundwater nitrate concentration , 3.0 mg L À1 as the indicator for assessing groundwater quality.Nitrate concentration above 3.0 mg L À1 is likely caused by anthropogenic activities, and this level is considered a contamination threshold (Tesoriero and Voss 1997).Specifically, we performed probability kriging analysis to interpolate nitrate data from 528 shallow groundwater wells (Appendix: Fig. D1) obtained from Groundwater Retrieve Network (GRN) and Wisconsin Department of Natural Resources (WDNR).We mapped the interpolation results at a 30-m spatial resolution using Geostatistical Analyst in ArcGIS 10.0 (Appendix: Fig. C1).Cross-validation was performed, and the result suggested that our kriged map predicted nitrate concentration across the watershed well (Qiu and Turner 2013; see Appendix).
where ALV x is the adjusted phosphorus export from pixel x, E y is the filtration efficiency of each downstream pixel y, and x represents phosphorus transport route from where it originated to the downstream water bodies, and ALV x was calculated as where pol x is the phosphorus export coefficient at pixel x, and HSS x is the hydrologic sensitivity score at pixel x, which was calculated as where logðR U Y U Þ calculated the sum of water yield of pixels along the flow path above pixel x (also called runoff index at pixel x), and d w is the calculated mean runoff index for the subwatershed of interest.We ran the model for 2006 at 30-m resolution for the watershed (Appendix: Fig. C1).The ecosystem service of providing high-quality surface water was the inverse of phosphorus loading-areas with lower phosphorus loading values delivered more surface-water quality service, and vice versa.We assessed our estimates by comparing with USGS gauge monitored phosphorus loading data for two subbasins of this watershed; our estimates agreed well with the gauge data (Qiu and Turner 2013; see Appendix).

Landscape pattern analysis
National Land Cover Data (NLCD; Fry et al. 2011) for 2006-the same year for which all hydrologic services were assessed-were used for landscape pattern analysis.NLCD data were derived from classification of Landsat Thematic Mapper imagery and have a spatial resolution of 30-m.Because aggregation improves the accuracy of individual land-cover classification and accuracy of calculated landscape metrics (Jones et al. 2001), we aggregated the 16 NLCD landcover data into seven broad classes (Table 3) prior to landscape pattern analysis.
Numerous landscape metrics have been developed to quantify composition and spatial configuration of cover types (Gustafson 1998, McGarigal et al. 2012).We selected a subset of metrics that are ecologically relevant (with hypotheses; Table 2) and may affect ecohydrological processes.Composition was quantified using the proportion of each land use/cover (PLAND); configuration metrics included: (1) fragmentation indices: patch density (PD) or edge density (ED); (2) connectivity indices: patch cohesion index (COHESION) and proportion of like adjacencies (PLADJ); (3) dispersion/interspersion indices: contagion index (CONTAG; Table 2).Most metrics were calculated by landcover type (class-level), except for contagion, which is calculated for the landscape.Please refer to McGarigal et al. (2012) for greater details on definition, computation and interpretation of each metric.All selected landscape metrics were calculated for each subwatershed using FRAG-STATS 4.0 (McGarigal et al. 2012).

Statistical analyses
Individual hydrologic ecosystem services were imported into ArcGIS 10.0 for visual representation and spatial pattern analysis.We evaluated the degree of spatial clustering (among subwatersheds) of each hydrologic service using Moran's I (with queen contiguity).
Because landscape composition constrains configuration (Gardner et al. 1987, Gustafson and Parker 1992, Gustafson 1998), we evaluated relationships between landscape pattern and hydrologic services [question 1] using a two-step procedure.The first step examined explanatory power of landscape composition variables alone for each service.As composition of all land  2).In addition, because percent cropland and urban were highly correlated (Pearson's r ¼ À0.9), we intentionally selected one or the other (never both) in the analysis, based on their hypothesized importance for the particular service.The second step used a priori models to test whether adding configuration variables for cover types whose composition was determined significant in the first-step analysis, could significantly improve the model after accounting for effect of composition.In both steps, we performed spatial regression for freshwater supply and surface-water quality because preliminary analyses indicated significant spatial autocorrelation in the residuals.We used multiple linear regression for groundwater quality, as no spatial structure was detected.For spatial regression, a spatial error model was chosen over a spatial lag model based on Lagrange multiplier tests (Bivand et al. 2008).We chose a spatial regression model with simultaneous autoregressive error term (i.e., SAR model) over a conditional autoregressive (CAR) model because SAR models had lower Akaike's information criteria (AIC) values and less remaining spatial autocorrelation in the residuals than equivalent CAR models (Bivand et al. 2008).To implement spatial models, we adopted first-order neighbor contiguity and computed spatial weights using row standardization.Rather than making inference from a single best model, we used multi-model inference analysis to avoid the uncertainty of model selection (Burnham and Anderson 2002).We used second-order AIC (i.e., AIC c ) corrected for small sample sizes to evaluate and rank all possible models, and created a subset of models within 2 AIC c units of the best model for model averaging (Burnham and Anderson 2002).The final output provided model-averaged coefficient estimates as well as the importance (i.e., the sum of the model weights within the set that included that variable) for each predictor variable.A Moran's I test for the full model indicated little spatial autocorrelation remained in the residuals after fitted with spatial regression (Table 4).We analyzed landscape effects on groundwater quality following the same procedure as described above but using linear regressions.We assessed the fit of full models using Nagelkerke pseudo R 2 for spatial regression and adjusted R 2 for linear regression, and compared between models with composition alone and their counterparts with composition and configuration variables, to determine the extent to which adding configuration variables can improve the model.All predictor variables were standardized prior to analyses.Heterogeneity of residuals, normality of errors, multicollinearity among predictor variables indicated by variance inflation factors (VIFs), and overly influential data points with Cook's distance were assessed for the full models, and no violations were detected.All statistical analyses were performed in the R statistical software (R Development Core Team 2007).We used OpenGeoDa (Anselin et al. 2005) and 'spautolm' function in 'spdep' package (Bivand 2010) for spatial regression, 'lm' function for multiple linear regression, and 'model.avg'function in the 'MuMIn' package (Barton 2010) for multi-model inference analysis.
To address our second research question, we first inspected scatterplots for indications of nonlinear responses in the provision of hydrologic services to landscape variables that were identified significant in the previous analyses.If a nonlinearity was suggested, we then tested for its significance by comparing models fitted with linear regression to models fitted using any of the three nonlinear functions (Fig. 1)-exponential curve, sigmoid logistic curve and step functionbased on AIC values (Burnham and Anderson 2002).Nonlinear responses were indicated if the nonlinear models had lower AIC values than linear models, and we then performed nonparametric deviance reduction analysis (Qian et al. 2003, King et al. 2005).Nonparametric deviance reduction analysis is designed to detect an ecological change-point along an environmental gradient (Qian et al. 2003, King et al. 2005).
Nonparametric deviance reduction analysis estimates the numerical value of a predictor x, resulting an abrupt change in the response variable, y, represented as the cumulative probability of a change-point.This method performs well if there is a sound rationale behind a revealed ecological change-point.Because this approach is based on deviance reduction, it will always detect a change-point along an environmental gradient regardless of whether there is a real ecological change.Hence, we performed the approximate v 2 test to determine whether the resulting ecological change-point was statistically significant, and also used a bootstrap simulation (N ¼ 1,000) to calculate 95% confidence interval (CI) of the estimated change-point.To illustrate the potential applicability of this approach, we identified subwatersheds in which minimal alterations (up to 5% of the watershed) in land cover would exceed the identified ecological change-point.We used 'nls' function in nls2 package for nonlinear regression, and the custom function 'ncpa' (Qian et al. 2003) for nonparametric deviance reduction analysis in the R statistical software (R Development Core Team 2007).

Spatial pattern of hydrologic ecosystem services among subwatersheds
Production of individual hydrologic services varied geographically among subwatersheds (Fig. 3; Appendix: Fig. C1).All three services were spatially aggregated rather than randomly distributed on the landscape (all Moran's I .0.31; P , 0.001).Estimates for freshwater supply, indicated by annual groundwater recharge, varied from 26.3 to 42.8 cm yr À1 among subwatersheds (Table 1).Groundwater quality, estimated by the probability of nitrate concentration ,3.0 mg L À1 , ranged from 0.05 to 0.92 among subwatersheds (Table 1).Estimates for surface-water quality, indicated by annual phosphorus loading, ranged from 0.07 to 0.17 kg ha À1 among subwatersheds (Table 1).

Landscape pattern effects on hydrologic ecosystem services
Freshwater supply.-Bothlandscape composition and configuration variables affected provision of freshwater supply, as indicated by the most supported models (Appendix: Table A1) and model-averaged parameter estimates (Fig. 4).Freshwater supply was negatively correlated with percent urban and wetland, and positively correlated with urban edge density (Fig. 4).There was a marginally significant relationship between freshwater supply and percent grassland (b ¼ 0.58; 95% CI ¼ À0.02-1.17,slightly overlapping zero; Fig. 4).Among all significant variables, percent wetland had the strongest effect on freshwater supply (b ¼ À1.72).Composition variables had a strong influence on freshwater supply and explained most of the variation (Nagelkerke pseudo R 2 ¼ 60%); configuration metrics explained an additional 4% of the variance (Table 4).
Groundwater quality.-Amongselected landscape variables, only percent cropland and grassland were statistically significant (all P , 0.001) for groundwater quality.Percent cropland and grassland were negatively correlated with groundwater quality (Fig. 4; Appendix: Table A2), and these two variables together explained 39% of the variation (Table 4).Comparatively, percent cropland had a stronger influence than grassland (b ¼ À0.10 vs. À0.06).The spatial configuration of cropland and grassland was not important for this service.

Surface-water quality
Landscape composition and configuration both influenced surface-water quality and explained more variation than for other two v www.esajournals.orgservices (Nagelkerke pseudo/adjusted R 2 ¼ 86%, vs. 64% for freshwater supply and 39% for groundwater quality).Suggested by the best models (Appendix: Table A3) and model-averaged coefficients (Fig. 4), greater surface-water quality (i.e., lower phosphorus loading) was associated with higher percentages of forest, grassland and wetland, and lower percentage of cropland.Spatial configuration of land covers also mattered.After accounting for composition, subwatersheds with higher wetland patch density, higher grassland patch density, more disaggregated forest patches and lower contagion had greater supply of surface-water quality service.The spatial arrangement of cropland (e.g., cropland edge density, cropland patch cohesion), v www.esajournals.orghowever, did not affect this indicator.Among all significant variables, percent cropland and forest had the most pronounced effects (b ¼ 0.013 and À0.008, respectively, for the indicator of phosphorus loading; Fig. 4).Landscape composition had a strong influence on surface-water quality and explained most of the variation (Nagelkerke pseudo R 2 ¼ 82%); adding class-and landscapelevel configuration variables explained an additional 4% of the variance (Table 4).

Nonlinear responses of hydrologic service to landscape pattern
We detected nonlinear relationships between surface-water quality and percent cropland and wetland (Fig. 5A, C).For both landscape variables, models fitted with an exponential function had the lowest AIC values (percent cropland: À91.5 vs. À84.4 with linear model; percent wetland: À49.1 vs. À46.5 with linear model).
Nonparametric deviance reduction analysis suggested that large changes in surface-water quality at the subwatershed scale occurred at 62% for percent cropland (95% CI: 60.3-66.3%;Fig. 5B), and 3% for percent wetland (95% CI: 1.9-5.8%;Fig. 5D).Given up to 5% land use/cover could be modified within a subwatershed, reducing cropland below 60% or increasing wetland above 6% might largely mitigate phosphorus loading and enhance surface-water quality.We also mapped subwatersheds that satisfied both criteria (Fig. 6) to identify locations where small land-cover change (converting 5% of the subwatershed area from cropland to wetland) could be prioritized to improve surface-water quality.

DISCUSSION
Results of this study provide empirical evidence for the role of landscape heterogeneity in v www.esajournals.orgsustaining three hydrologic ecosystem services in an urbanizing agricultural watershed.Understanding effects of landscape pattern on ecosystem services is an important research goal and a basis from which tradeoffs, synergies, trajectories and management alternatives can be considered (Turner et al. 2013).We found a consistently stronger influence of landscape composition than of spatial configuration, and our analyses suggested potential thresholds for land-cover proportions related to surface-water quality.While some of the relationships have been observed or modeled in earlier studies (Soranno et al. 1996, Jones et al. 2001, Strayer et al. 2003), this research expands understanding of landscape pattern effects by simultaneously considering multiple hydrologic services and explicitly quantifying effects of configuration.

Does spatial pattern matter?
Our results suggested that landscape composition and configuration both mattered for supply of hydrologic services, but composition had by far the stronger influence (Fig. 4, Table 4).Other empirical or modeling studies also point to a dominant influence of composition over spatial configuration.For example, in the literature on effects of habitat fragmentation on biodiversity, habitat loss (i.e., composition) shows a consistently stronger effect than habitat fragmentation per se (i.e., configuration;Fahrig 2002Fahrig , 2003)).A global synthesis of pollination also revealed that the amount of habitat was more important than spatial configuration (Kennedy et al. 2013).Theoretical studies suggest that landscape composition in general exerts a stronger influence, and the effect of landscape configuration should be most apparent when land-cover types are present at intermediate levels of abundance (Gergel 2005).When the landscape is dominated (e.g., .60%)by a given cover type, it is nearly always well connected; when a cover type is scarce, it is never well connected (Gergel 2005).Our results provide empirical support for this argument and suggest where agricultural or urban cover dominate subwatersheds, landscape Fig. 6.Spatial distribution of subwatersheds (indicated in red) where relatively small changes in land cover (i.e., converting 5% of the area from cropland to wetland) could largely enhance surface-water quality.
v www.esajournals.orgcomposition has a greater effect on hydrologic services and configuration matters less.
Groundwater recharge, indicative of freshwater supply, was negatively correlated with percent urban and wetland in subwatersheds.Urban cover, primarily consisting of impervious surface, increases runoff and reduces infiltration that replenishes groundwater (Arnold and Gibbons 1996).However, the negative effect of wetland on recharge was surprising.A synthesis suggests that wetlands can enhance or impede recharge, depending on geographic location and soils of wetland (Bullock and Acreman 2003).In our study area, many wetlands are water sources (i.e., springs) occupying low topographic position where groundwater is discharged (Winter 1999).Other wetlands overlie impermeable soils or rocks that hinder infiltration and groundwater recharges.Our results also suggested that for a given amount of urban area, increased urban edge density enhanced groundwater recharge.Increasing the interface of urban and more vegetated covers, as occurs with urban green space, may reduce the rate of surface flow and promote interactions of water with soils and vegetation, thus increasing infiltration of water.
The effect of percent cropland on groundwater quality agrees with earlier studies identifying intensive agriculture as the main cause of elevated groundwater nitrate levels (Spalding and Exner 1993).In this watershed where corn and dairy farms are dominant, nitrogen inputs from manure and fertilizer largely exceed plant demand, leading to increased nitrate leaching.The negative relationship of percent grassland and groundwater quality is probably because our grassland category included forage crops (Fig. 2) that receive nitrogen fertilizer, albeit in smaller amounts.Further, forage crops and row crops can be difficult to distinguish from each other in remote sensing imagery classification.It is also possible that grassland in 2006 may have been transitioned from cropland in previous years.Thus, misclassification or lag effects might have contributed to the detected relationship between percent grassland and groundwater quality.The lack of relationship between configuration variables and groundwater quality is consistent with the fact that nitrate is very soluble and can move quickly along vertical flows.Overall, landscape variables explained only 39% of the variation (Table 4), suggesting other factors such as soil, geology and climate may play an important role in affecting groundwater quality.
Surface-water quality was negatively correlated with percent cropland and positively correlated with percent forest, wetland and grassland.This result is consistent with previous findings that fertilizer and manure in agriculture are major sources of phosphorus loading (Carpenter et al. 1998).In contrast, natural vegetation, especially as riparian buffers, can attenuate sediment and reduce nutrient flows (Peterjohn andCorrell 1984, Soranno et al. 1996).Landscape configuration also influenced phosphorus loads, likely because phosphorus adsorbs to sediment, and lateral fluxes of sediment and runoff can be mediated by land-cover arrangement.The effect of wetland patch density suggests that many small and widely distributed wetlands increase opportunities to trap nutrients (Knight 1992).Similarly, forests appeared to be more effective for retaining nutrients when dispersed on the landscape.At the landscape level, lower contagion indicates fine-scale interspersion of land covers, which may facilitate nutrient attenuation.Spatial configuration of cropland, however, was unimportant, perhaps because cropland is so dominant in many of these subwatersheds.
When comparing among services, landscape pattern had varying or even opposing effects on different hydrologic services (Fig. 4), resulting in tradeoffs.For example, wetland restoration could improve surface-water quality, but could reduce recharge.Grassland restoration might enhance surface-water quality and groundwater recharge (though marginally), but could lead to deterioration of groundwater quality if grasslands were fertilized.Prior research has reported similar empirical patterns of ecosystem service tradeoffs or synergies (Raudsepp-Hearne et al. 2010, Qiu andTurner 2013).Multiple consequences of landscape pattern for hydrologic services suggest that improving the provision of all services is likely to be challenging.

Small change, large gains?
Our study suggested opportunities to enhance surface-water quality by reducing percent cropland below 60% and/or increasing wetland above 6% in subwatersheds.Percolation theory may explain why it is important to reduce cropland to v www.esajournals.orgbelow 60% (Turner et al. 2001, Gergel 2005).When the proportion of cropland exceeds ;60%, croplands coalesce into larger and more connected patches, therefore reducing edge effects and nutrient retention by other land covers.Jones et al. (2001) implied that the amount of agriculture may reach a threshold for nutrient loads, and our findings support their hypothesis.While the rationale for the wetland threshold remains less clear, nonetheless, it suggests the functional importance of wetlands to exceed this relative abundance in a subwatershed.
Nonlinear responses of surface-water quality to percent cropland and wetland suggest that, given limited amount of allowable land-use/ cover change, choosing subwatersheds to change does matter.For example, reducing cropland from 65% to 55% in a subwatershed may have a greater effect on enhancing surface-water quality than reducing cropland from 45% to 35% (Fig. 5).Similarly, restoring wetland from 1% to 6% may provide greater benefits for surface-water quality than restoring wetland from 15% to 20% (Fig. 5).Although these values are likely context dependent and should be evaluated in other landscapes, our study suggests that understanding nonlinear responses of ecosystem services to landscape pattern could help identify areas where small management may benefit ecosystem service supply.

Management implications
Demand for scientific rationales to manage landscapes is substantial (Turner 2005).Understanding how landscape patterns influence ecosystem services will have important implications for planning, managing and restoring landscapes.We found that in agricultural-dominated landscapes, composition had prevailing influences relative to spatial configuration, suggesting that altering spatial arrangement of land cover alone will be inadequate to improve hydrologic services within subwatersheds.Rather, the relative amount of land-cover classes may need to change to improve hydrologic services.Spatial configuration can mediate services dependent on lateral transfers, and our analysis showed the importance of edge effects (e.g., urban or forest edges) for freshwater supply and surface-water quality.Where land-cover change is not feasible (e.g., highly urbanized areas), there might be opportunities to create more edges with vegetated cover types such as parks, conservancy lands and rain gardens, to enhance hydrologic services (Brosi et al. 2008).In addition, the magnitude of effects varied among land-cover classes, suggesting the importance of targeting cover types with stronger effects for management.Cropland, for instance, had greater effects on water-quality services than other cover types (Fig. 4); thus reducing croplands and associated nutrient inputs may provide more benefits in mitigating nitrate leaching and phosphorus loads to freshwaters.Overall, landscape pattern will be a more effective tool for enhancing surface-water quality than for groundwater quality or freshwater supply.
Understanding nonlinear responses has important management implications, in particular if the extent to which landscapes can be altered is limited.Conservation planning through land easement, acquisition and incentives, for example, is constrained by resource availability, and the amount of land that can be conserved is often very limited.Hence, priorities should be placed on locales where per unit land-use/cover change on enhancing ecosystem services is greater.By strategically targeting a small number of subwatersheds in this watershed (Fig. 6), we may improve hydrologic services for least cost.However, more research is needed to identify other variables that can also lead to nonlinear responses, and expand to other services that are also essential for human well-being.

Methodological considerations
We predicted indicators of freshwater supply and surface-water quality in 30 3 30-m grid cells using spatially explicit simulation models in which land cover was among the inputs.Although land-cover proportions are subsequently included at the sub-watershed level to test for effects of landscape composition and configuration, there is no a priori deterministic relationship between land cover and the indicators (groundwater recharge and phosphorus loading) we modeled.Land cover interacts with soil, climate and topography to predict indicators in each grid cell, and it is only involved in a subset of the simulated biophysical processes.Second, the scales at which the indicators are simulated and landscape patterns are analyzed differ, and we v www.esajournals.orgaimed to determine whether effects of land cover were detectable at the sub-watershed scale after grid-cell processes were simulated.
Our study focused on disentangling effects of land cover composition and configuration in a region with minimal topographic relief, and we did not consider topographic or hydrologic gradients.Our study region was glaciated, and the terrain is generally flat (mean slope of 2.28, with standard error of 0.128 for the subwatersheds).Topographic and hydrologic gradients may be of much greater importance in areas with mountains or complex topographies.However, hydrologic gradients are likely to influence groundwater-related indicators, and these warrant future research.

CONCLUSIONS
Our study revealed a consistently dominant influence of composition relative to spatial configuration for hydrologic services, and suggested potential nonlinear responses of surfacewater quality to landscape pattern.Relationships observed in this study are applicable to other agriculture-dominated watersheds with little topographic relief, but results could differ in landscapes with less human influence, or where topographic or other environmental gradients are much stronger.Our study contributed to an emerging literature on understanding linkages between landscape pattern and the provision of ecosystem services (Turner et al. 2013).Knowledge gained from this study has the potential to enhance management of agricultural landscapes for sustaining the supply of hydrologic ecosystem services and creating multifunctional landscapes.

Fig. 1 .
Fig. 1.Hypothetical nonlinear responses of ecosystem service to landscape pattern: (A) logistic or exponential curve; (B) sigmoid logistic curve; (C) step function.Solid purple circles indicate areas approaching the threshold and suggest where small changes in landscape pattern can yield substantial gains in ecosystem services provision.

Fig. 2 .
Fig. 2. Map of the Yahara River Watershed (Wisconsin, USA) and the land use/land cover pattern (with percent cover) for 2006, derived from National Land Cover Data (NLCD).Land use/cover was aggregated into seven broad categories for landscape pattern analysis.Delineations of the Yahara Watershed and subwatersheds (n ¼ 100) were based on stream network data, light detection and ranging (LiDAR) elevation, and a field-checked basin map.Two aerial photos were taken in summer 2013 and illustrated typical agricultural landscape (upper) and increasing urbanization (lower) in this watershed.

Fig. 3 .
Fig. 3. Spatial distributions of three hydrologic ecosystem services: (A) freshwater supply; (B) groundwater quality, and (C) surface-water quality, at the subwatershed level for 2006.Green indicates areas with high supply and purple indicates low supply of hydrologic services.

Fig. 4 .
Fig. 4. Model-averaged coefficient estimates (b) for effects of landscape composition (orange) and configuration (green) on the biophysical indicators (in parentheses) of all three hydrologic ecosystem services.The bar around the mean indicates 95% confidence intervals (CI).The landscape pattern effect is significantly different from zero when its 95% CI do not bracket zero.

Table 2 .
Landscape composition and configuration variables and their hypothetical relationships (Rel; ''þ'' denotes positive and ''À'' denotes negative relationships) with the biophysical indicators (in parentheses) of each hydrologic service.

Table 3 .
Land use/land cover reclassification used in this study compared with original NLCD 2006 landcover classes.

Table 4 .
Summary of the full models for landscape pattern effects on three hydrologic services: variance explained by landscape composition alone vs. explained by landscape composition and configuration ( R 2 denotes Nagelkerke R 2 ); Moran's I in model residuals; and variance inflation factors (VIF) for testing multicollinearity among predictor variables.