Methods

Each ecosystem service was quantified and mapped by using empirical estimates and spatially explicit model for the terrestrial landscape of the Yahara Watershed for 2006. Crop production (expected annual crop yield, bu per yr) Crop yield was estimated for the four major crop types (corn, soybean, winter wheat and oats) that account for 98.5 percent of the cultivated land in the watershed by overlaying maps of crop types and soil-specific crop yield estimates. The spatial distribution of each crop was obtained from the 2006 Cropland Data Layer (CDL) from the National Agricultural Statistics Service (NASS) and soil productivity data were extracted from Soil Survey Geographic (SSURGO) database. Crop and soil data were converted to 30 m resolution and the two maps were overlain to estimate crop yield in each cell. For each crop-soil combination, crop area was multiplied by the estimated yield per unit area. Estimates for each crop type were summed to map estimated crop yield for 2006. Pasture production (expected annual forage yield, animal-unit-month per year ) As for crop production, forage yield was estimated by overlaying the distribution of all forage crops (alfalfa, hay and pasture/grass) and soil specific yield estimates. The spatial distribution of each forage crop was also derived from 2006 CDL, and rescaled to 30 m grid prior to calculation. The SSURGO soil productivity layer provided estimates of potential annual yield per unit area for each forage crop. Overlay analyses were performed for each forage-soil combination, as done for crops, and summed to obtain the total expected forage yield in the watershed for 2006. Freshwater supply (annual groundwater recharge, cm per year) . Groundwater recharge was quantified and mapped using the modified Thornthwaite-Mather Soil-Water-Balance (SWB) model. SWB is a deterministic, physically based and quasi three-dimensional model that accounts for precipitation, evaporation, interception, surface runoff, soil moisture storage and snowmelt. Groundwater recharge was calculated on a grid cell basis at a daily step with the following mass balance equation<p align="center">Recharge= (precipitation + snowmelt + inflow) &ndash;<p align="center">(interception + outflow + evapotranspiration) &ndash; delta soil moisture<p align="center"> We ran the model for three years (2004 to 2006) at 30m resolution, with the first two years as spin up of antecedent conditions (e.g. soil moisture and snow cover) that influence groundwater recharge for the focal year of 2006. Carbon storage (metric tonsper ha) We estimated the amount of carbon stored in each 30 m cell in the Yahara Watershed by summing four major carbon pools: aboveground biomass, belowground biomass, soil carbon and deadwood/litter. Our quantification for each pool was based mainly on carbon estimates from the IPCC tier-I approach and other published field studies of carbon density and was estimated by land-use/cover type.Groundwater quality (probability of groundwater nitrate concentration greater than 3.0 mg per liter, unitless 0 to1) Groundwater nitrate data were obtained from Groundwater Retrieve Network (GRN), Wisconsin Department of Natural Resources (DNR). A total of 528 shallow groundwater well (well depth less than the depth from surface to Eau Claire shale) nitrate samples collected in 2006 were used for our study. We performed kriging analysis to interpolate the spatial distribution of the probability of groundwater nitrate concentration greater than 3 mgper liter. We mapped the interpolation results at a 30m spatial resolution using Geostatistical Analyst extension in ArcGIS (ESRI). In this map, areas with lower probability values provided more groundwater quality service, and vice versa. Surface water quality (annual phosphorus loading, kg per hectare). We adapted a spatially explicit, scenario-driven modeling tool, Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) to simulate discharge of nonpoint-source phosphorus. A grid cells phosphorus contribution was quantified as a function of water yield index, land use/cover, export coefficient, and downslope retention ability with the following equation:Expx = ALVx * sum of the products from y=x+1 to X for (1-Ey)where ALVx is the adjusted phosphorus export from pixel x , Ey is the filtration efficiency of each downstream pixel y , and X represents phosphorus transport route from where it originated to the downstream water bodies. Filtration efficiency was assigned by cover type: natural vegetation was assigned a high value, semi-natural vegetation an intermediate value, and developed or impervious covers were assigned low values. We ran the model for 2006 and mapped estimated phosphorus loading across the watershed. The ecosystem service of providing high quality surface water was the inverse of phosphorus loading. Therefore, areas with lower phosphorus loading values delivered more surface water quality, and areas with higher phosphorus loading values supplied less surface water quality.Soil retention (annual sediment yield, metric tons per hectare). We quantified annual sediment yield as the (inverse) indicator for soil retention by using the Modified Universal Soil Loss Equation (MUSLE). MUSLE is a storm event based model that estimates sediment yield as a function of runoff factor, soil erodibility, geomorphology, land use/cover and land management. Specifically, a grid cells contribution of sediment for a given storm event is calculated as:Sed= 11.8*(Q*q_p)^0.56 * K * LS * C * Pwhere Sed represents the amount of sediment that is transported downstream network (metric tons), Q is the surface runoff volume (m³), q_pis the peak flow rate (cubic meters per s), K is soil erodibility which is based on organic matter content, soil texture, permeability and profiles, LS is combined slope and steepness factor, and C* P is the product of plant cover and its associated management practice factor. We used the ArcSWAT interface of the Soil and Water Assessment Tool (SWAT) to perform all the simulations. We ran this model at a daily time step from 2004 to 2006, with the first two years as spin up , then mapped total sediment yield for 2006 across the watershed. Similar to surface water quality, the ecosystem service of soil retention was the inverse of sediment yield. In this map, areas with lower sediment yield provided more of this service, and areas with higher sediment yield delivered less. Flood regulation (flooding regulation capacity, unitless, 0 to 100) We used the capacity assessment approach to quantify the flood regulation service based on four hydrological parameters: interception, infiltration, surface runoff and peak flow. We first applied the Kinematic Runoff and Erosion (KINEROS) model to derive estimates of three parameters (infiltration, surface runoff and peak flow) for six sampled sub basins in this watershed. KINEROS is an event-oriented, physically based, distribution model that simulates interception, infiltration, surface runoff and erosion at sub-basin scales. In each simulation, a sub basin was first divided into smaller hydrological units. For the given pre-defined storm event, the model then calculated the amount of infiltration, surface runoff and peak flow for each unit. Second, we classified these estimates into 10 discrete capacity classes with range from 0 to 10 (0 indicates no capacity and 10 indicates the highest capacity) and united units with the same capacity values and overlaid with land cover map. Third, we calculated the distribution of all land use/cover classes within every spatial unit (with a particular capacity). We then assigned each land use/cover a capacity parameter based on its dominance (in percentage) within all capacity classes. As a result, every land use/cover was assigned a 0 to 10 capacity value for infiltration, surface runoff and peak flow. This procedure was repeated for six sub basins, and derived capacity values were averaged by cover type. We applied the same procedure to soil data and derived averaged capacity values for each soil type with the same set of three parameters. In addition, we obtained interceptions from published studies for each land use/cover and standardized to the same 0 to 10 range. Finally, the flood regulation capacity (FRC) for each 30m cell was calculated with the equation below:FRC= for each land use and land cover class the sum of (interception + infilitration + runoff + peakflow) + for each soil class the sum of (infiltration + runoff + peakflow).To simplify interpretation, we rescaled original flood regulation capacity values to a range of 0 to100, with 0 representing the lowest regulation capacity and 100 the highest. Forest recreation (recreation score, unitless, 0 to 100). We quantified the forest recreation service as a function of the amount of forest habitat, recreational opportunities provided, proximity to population center, and accessibility of the area for each 30m grid cell with the equation below:FRSi= Ai * sum of (Oppti + Popi + Roadi)where FRS is forest recreation score, A is the area of forest habitat, Oppt represents the recreation opportunities, Pop is the proximity to population centers, and Road stands for the distance to major roads. To simplify interpretation, we rescaled the original forest recreation score (ranging from 0 to 5200) to a range of 0 to 100, with 0 representing no forest recreation service and 100 representing highest service. Several assumptions were made for this assessment approach. Larger areas and places with more recreational opportunities would provide more recreational service, areas near large population centers would be visited and used more than remote areas, and proximity to major roads would increase access and thus recreational use of an area. Hunting recreation (recreation score, unitless 0 to100) We applied the same procedure used for forest recreation to quantify hunting service. Due to limited access to information regarding private land used for hunting, we only included public lands, mainly state parks, for this assessment. The hunting recreation service was estimated as a function of the extent of wildlife areas open for hunting, the number of game species, proximity to population center, and accessibility for each 30m grid cell with the following equation:<br />HRSi= Ai * sum of (Spei + Popi + Roadi)where HRS is hunting recreation score, A is the area of public wild areas open for hunting and fishing, Spe represents the number of game species, Pop stands for the proximity to population centers, and Road is the distance to major roads. To simplify interpretation, we rescaled the original hunting recreation score (ranging from 0 to 28000) to a range of 0 to100, with 0 representing no hunting recreation service and 100 representing highest service. Similar assumptions were made for this assessment. Larger areas and places with more game species would support more hunting, and areas closer to large population centers would be used more than remote areas. Finally, proximity to major roads would increase access and use of an area.

Methods

SamplingAt each location, the extent of the peatland basin was examined visually using the soils map and a long axis was defined as the longest linear stretch of peat, while a short axis was defined perpendicular to the long axis. The sample area of a given site was defined as the entire peatland basin (full basin site, N = 11) if the length of the long axis was 800 m or less. For larger peatlands, the site was defined as an area with width 150-200 m and length 400-600 m (partial basin site, N = 10) extending outward from one edge of the peatland and if possible crossing the entire short axis. Peat depth was measured throughout the area on a regular grid at intervals varying from 20-90 m depending upon the size of the site. In addition, vegetation was surveyed and peat pore water chemistry was sampled at 3 plots located at 25%, 50% and 75% of the length of the long axis of the sampling area. Peat cores were taken at the same plots for a subset of 5 sites described below, and slope at the upland-peatland interface at the edge of the site was also measured for all 11 full basin sites and 4 of the 10 partial basin sites.The depth of organic sediment (primarily peat) was measured to depth of contact with mineral surface (typically sand) throughout the sampling area using a stainless steel peat depth probe (PDP) on a regular grid at intervals varying from 20-90 m. Two different versions of the PDP were used, and intercalibrated to ensure consistency. The first consisted of 60 (1.83 m) sections of 3 800 (0.95 cm) diameter threaded steel rod, connected with hex-shaped coupling nuts. The second was a custom-made version with the same general design including length and diameter of sections, but consisted of a smooth stainless steel surface and contained an inset male and female threading system to avoid the protruding coupling nuts. The PDP was used only to determine depth to refusal and was not equipped to collect samples; thus it could not differentiate between peat and soft organic lacustrine sediment. In nearly all cases, the person using the PDP could feel contact with sand (typical glacial sediments) at depth to refusal.Peat CoresPeat cores were taken at 13 different locations, including the central plot for site 4n and each of the 3 plots for sites 7b, 9b, 12f, and 21b. At each core location, samples were taken using a Russian-style corer (50 cm length, 5cm diameter) at depths of 0.5, 1.0, 2.0, 4.0 and 6.0 m, up to the maximum depth (peat-sand interface). We examined peat color and degree of decomposition using the von Post scale in the field. Particular attention was paid to the presence/absence of gyttja at the peat-sand interface. Gyttja is a dark olive-green algae-derived gelatinous lacustrine sediment, which indicates the former presence of a clear-water lake at a given site.For each core sample (N = 45), a central 10 cm section was preserved and used to measure moisture content, bulk density, and organic matter (OM) content in the laboratory. The 10 cm section was halved vertically, and one half (between 50 and 150 g wet weight) was used for measurement of wet bulk density (rw =mw per V; where mw = wet weight and V = volume measured by water displacement). The other half was used to measure mass loss by oven-drying at 55C until the mass was stable (typically 5-10 days, measurement precision = or - 0.1 g). Volumetric moisture content was calculated as (mw md) per V and bulk density (rb)asmdper V where md =dry weight, mw = wet weight, and V = volume calculated as mw per rw. From the dried sample, a 1-3 g homogenized subsample was ashed in a muffle furnace at 440C for 8 h to determine ash-free dry weight (maf = md mash; precision plus or minus 0.01 g), and OM content (OM%) was calculated as maf per md. Finally, OM density (rOM) was calculated as rb OM%.For each of the 13 core locations, we estimated the total mass of OM by summing the product of rOM and volume over all measurement intervals. To estimate a continuous vertical distribution, rOM was interpolated linearly by depth between measurement points. The 0.25-0.5 m interval was assigned the same rOM as the 0.5 m value, while the 0-0.25 m surface interval was assigned a rOM of half of that measured at 0.5 m, to account for the lower bulk density in living recently dead Sphagnum in the acrotelm. The deepest measured rOM value was extrapolated down to a depth of 0.25 m above the base of the core, and the basal 25 cm of the core was assigned a rOM of 46 kg m3. This is equivalent to the mean value measured for gyttja, to account for the fact that the peat is grading into lower-OM gyttja and or sand at the interface with glacial till. Vertically averaged mean rOM was calculated as the total mass of OM in the core divided by the total core volume.Edge Slope in the FieldBecause many peatlands in this region formed from in-filling of lakes, we hypothesized that local geomorphology, specifically slope at the peatland margin (peatland-upland interface), might be a good indicator of peatland depth. At a subset of 15 sites (including all 11 full basin sites) we measured slope at the peatland-upland interface (Edge Slope in the Field, ESF). At full basin sites, ESF readings were taken at 8 peatland-upland interface locations distributed evenly around the edge of the site. At partial basin sites, measurements were only taken at those site edges that were adjacent to upland, resulting in fewer than 8 locations at each site. At each location, a Suunto clinometer was used to measure slope (%) from the peatland-upland interface oriented up the steepest upland slope at a distance of 5 m, 10 m, 20 m and 30 m, and these four values were averaged to give a single slope value for each location. The precision of individual measurements was plus or minus 0.3% slope (mean SD of replicate measurements). The values used for statistical analysis were site mean (ESFmean) and maximum (ESFmax) of location slopes.VegetationWe hypothesized that the surface vegetation characteristics might be related to peat depth, either directly (due to differential contributions of plant species to decomposition rates and water-holding capacity), or indirectly by responding to local environmental characteristics (e.g., water table, groundwater flow) that also influence peat formation. Vegetation was surveyed following a modification of the U.S. Forest Service Forest Inventory and Analysis (FIA) protocol. At each plot, a circular sampling area with 7.3 m radius was laid out, with 3 linear transects extending from the center to the perimeter at 0, 120 and 240 degrees from compass north. Within the circular plot, all trees with diameter at breast height (DBH) of 2.5 to 4.9 cm were counted as saplings and their species recorded, while species and DBH were recorded for all trees with DBH at 5 cm. Basal area for each tree species in units of m2 ha-1 was calculated by summing the DBH of all individual trees and normalizing by plot area. The mean height and intersection length of each of 4 categories of shrub (alder, bog birch, ericaceous or tree seedling) were recorded for woody vegetation of height greater than 50 cm (but DBH less than 2.5 cm) that intersected any of the linear transects. Shrub percent cover was estimated for each category by dividing the intercepted length by the total transect length. Coarse woody debris (CWD) with length greater than 1 m that intersected any of the linear transects with diameter greater than 5 cm was tallied; small end diameter (down to 5 cm), large end diameter, and length of each piece of CWD was recorded and used to calculate volume as described by Waddell [2002]. Finally, three ground-layer quadrats (1 m2) were laid out, 1 each adjacent to the 3 linear transects spanning a distance of 4 to 5 m from the center point. Within each quadrat, percent cover was recorded for each of 8 commonly occurring ground cover types: bare ground, ericaceous shrubs, ferns, forbs, graminoids, Sphagnum mosses, other mosses, and other woody vegetation (tree seedlings). Values for the 3 quadrats were averaged to give a plot mean cover of each ground cover type.

Methods

Riparian buffer strips within the Lake Mendota watershed, created by buffering the Lake Mendota feature in the Wisconsin Department of Natural Resources waterbody layer.