A coupled terrestrial and aquatic biogeophysical model of the Upper Merrimack River watershed , New Hampshire , to inform ecosystem services evaluation and management under climate and landcover change

Accurate quantification of ecosystem services (ES) at regional scales is increasingly important for making informed decisions in the face of environmental change. We linked terrestrial and aquatic ecosystem process models to simulate the spatial and temporal distribution of hydrological and water quality characteristics related to ecosystem services. The linked model integrates two existing models (a forest ecosystem model and a river network model) to establish consistent responses to changing drivers across climate, terrestrial, and aquatic domains. The linked model is spatially distributed, accounts for terrestrial–aquatic and upstream– downstream linkages, and operates on a daily time-step, all characteristics needed to understand regional responses. The model was applied to the diverse landscapes of the Upper Merrimack River watershed, New Hampshire, USA. Potential changes in future environmental functions were evaluated using statistically downscaled global climate model simulations (both a high and low emission scenario) coupled with scenarios of changing land cover (centralized vs. dispersed land development) for the time period of 1980–2099. Projections of climate, land cover, and water quality were translated into a suite of environmental indicators that represent conditions relevant to important ecosystem services and were designed to be readily understood by the public. Model projections show that climate will have a greater influence on future aquatic ecosystem services (flooding, drinking water, fish habitat, and nitrogen export) than plausible changes in land cover. Minimal changes in aquatic environmental indicators are predicted through 2050, after which the high emissions scenarios show intensifying impacts. The spatially distributed modeling approach indicates that heavily populated portions of the watershed will show the strongest responses. Management of land cover could attenuate some of the changes associated with climate change and should be considered in future planning for the region.


INTRODUCTION
Understanding changes in ecosystem services at regional scales is increasingly important for environmental scientists, managers, planners, and decision makers as climate and land-use change continue (de Groot et al. 2010).Management decisions require an improved understanding of the drivers and processes that influence ecosystem services.Change in major environmental drivers, such as climate and land cover, typically result in large changes in ecosystem service supply (Schroeter et al. 2005).These changes are interactive and complex across space and time (Chen et al. 2013), requiring the development of appropriate methods to elucidate functional tradeoffs between management strategies.In order for stakeholders and citizens to be able to assess the value of the ecosystem services being provided, they need them to be expressed in terms of indicators that clearly relate to environmental condition (Nelson et al. 2009, Carpenter et al. 2015, Qiu and Turner 2015).Indicators of ecosystem services should be quantifiable, scalable (Bagstad et al. 2013a, Carpenter et al. 2015), explicit in time and space, and sensitive to land-cover or management change (Burkhard et al. 2012, van Oudenhoven et al. 2012).Moreover, the appropriate indicator is dependent on the method by which the ecosystem service is being valued (de Groot et al. 2010).Commonly used indicators that are easily monetizable may be incomplete for more comprehensive sociocultural preference valuations (Mavrommati et al. 2017).
Aquatic ecosystems are strongly influenced by terrestrial environments in their watersheds.Indicators of environmental condition related to climate and land cover can be generated at regional scales from existing meteorological observations, downscaled projections, land-cover atlases, or land-cover change scenarios (Queiroz et al. 2015).However, these simple evaluations are difficult in freshwater systems because responses depend on terrestrial and climate conditions integrated over entire watersheds and that vary over time.Terrestrial environments, interacting with climate, determine flow, temperature, and nutrient regimes in regional drainage networks (Poff et al. 1997).Hydrology is the overriding control of many aquatic ecosystem services (Brauman et al. 2007) influencing water availability, instream habitats, water temperatures, and nutrient fluxes.Aquatic ecosystems are further influenced by internal processes such as temperature re-equilibration and nutrient removal as water flows from upstream to downstream (Hale et al. 2014).To understand changes in ecosystem services, linked terrestrial and aquatic ecosystem models are needed to capture the processes defining responses to changing land cover and climate, partition Ecology and Society 22(4): 18 https://www.ecologyandsociety.org/vol22/iss4/art18/ the influence of each, and permit estimation of responses beyond previously observed ranges.
A variety of spatially explicit models or tools have been used for watershed-scale studies of environmental indicators and ecosystem services.Examples include Artificial Intelligence for Ecosystem Services (ARIES) (Villa et al. 2009), Multiscale Integrated Models of Ecosystem Services (MIMES) (Boumans et al. 2015), and Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) (Tallis and Polasky 2009, Bagstad et al. 2013b, Tallis et al. 2013).These models do not capture seasonal or subseasonal climate variability, which is projected to change regionally (Wood et al. 2002, Hayhoe et al. 2007, Horton et al. 2014) and is important for capturing watershed functions related to flood attenuation, water provisioning, river temperature regulation, and other ecosystem services (Vigerstol and Aukema 2011).To fully account for changes in ecosystem function associated with altered precipitation, temperature, and land-use and land-cover patterns, process-based models that incorporate key space-and time-varying hydrological and ecological processes are critical (Bagstad et al. 2013b).
Several processed-based terrestrial and/or aquatic biogeophysical models have recently been used for ecosystem service valuation (Logsdon and Chaubey 2012a, Bagstad et al. 2013a, Carpenter et al. 2015).The Variable Infiltration Capacity (VIC) model is a large-scale, semidistributed hydrological model (Liang et al. 1994) that simulated provisioning hydrological ecosystem services (Vigerstol and Aukema 2011) and flood regulation (Lee et al. 2015).The Soil Water Assessment Tool (SWAT), a process-based, spatially distributed hydrological and water quality model (Notter et al. 2012) was used to evaluate aquatic environmental variables, including water yield (Karabulut et al. 2015) and water quality (Logsdon and Chaubey 2012b).To compute dynamic ecosystem services in the agricultural Yahara watershed, a linked terrestrialaquatic model used a process-based agroecosystem model (Agro-IBIS) (Soylu et al. 2014, Carpenter et al. 2015), a terrestrial hydrology model (THMB) (Coe 2000), a three-dimensional groundwater flow model (MODFLOW) (Harbaugh 2005), and a hydrological routing model (HYDRA), which lacked instream biogeochemistry (Coe 2000).There are a few robust examples of coupled human and biogeophysical models that have quantified ecosystem services at global (Boumans et al. 2002) and watershed (Costanza et al. 2002) scales.The Global Unified Metamodel of the Biosphere (GUMBO) simplifies several existing dynamic global models of both natural and social systems at an intermediate level of complexity and annual time-step (Boumans et al. 2002).The Patuxent Landscape Model (PLM; Costanza et al. 2002) is a spatially explicit process-based model that addresses the effects of both the magnitude and spatial patterns of human settlements and agricultural practices on hydrology, plant productivity, and nutrient cycling in the landscape, also at an annual time-step.Finally, valuation of ecosystem service information for land-use decisions has been estimated directly through an agent-based modeling framework (Heckbert et al. 2010, 2014, Groeneveld et al. 2017).
Our study contributes to the evolution of ecosystem service models by emphasizing shorter-term temporal dynamics in a spatially distributed and process-based framework that links terrestrial and aquatic function.This framework provides a new perspective for understanding the impacts that climate and landcover and land-use change have on terrestrial and aquatic resources.
We present an approach that links time-varying (daily time-step) terrestrial and aquatic ecosystem models at regional scales and apply this model into the future using scenarios of climate and land cover to project changes in ecosystem services.First, we describe an indicator framework that succinctly represents a comprehensive suite of environmental conditions relevant to important ecosystem services.Second, we describe the linkage and validation of the terrestrial and aquatic ecosystem models to simulate aquatic indicators through the 21st century.We integrated the Photosynthetic Evapotranspiration-Carbon and Nitrogen (PnET-CN) forest ecosystem model (Ollinger et al. 2002, 2008, Aber et al. 2005) and the Framework for the Aquatic Modeling of the Earth System (FrAMES) aquatic ecosystem model (Wollheim et al. 2008a, b, Wisser et al. 2010, Stewart et al. 2011, 2013, Mineau et al. 2015;Zuidema, Wollheim, Mineau, et al., unpublished manuscript).These models integrate the dynamics of terrestrial and aquatic processes and linkages at daily timesteps, making them ideal for studying aquatic ecosystem responses in forest-dominated watersheds.In coordination with a separate effort described elsewhere in this special issue (Mavrommati et al. 2017) to assess the value of ecosystem services provided by the Upper Merrimack River watershed (UMRW) of New Hampshire, we contrast two extremes of projected futures in climate and land-cover change.The outcome suggests that climate change influences most indicators of environmental condition in the UMRW more than changes in land cover, although land cover has important interactive capacity to dampen or exacerbate the effects of the changing supply of ecosystem services in the future.

METHODS
The central goal of this study was to develop estimates of ecosystem function critical to current and future watershed residents of the UMRW.We present a description and rationale for studying this watershed and the development of environmental indicators relevant to ecosystem services of local residents.We next describe the integrated terrestrial-aquatic model that projects aquatic environmental indicators, including scenarios, parameterizations, and data sets needed to project aquatic indicators to 2100.Finally, to understand how aquatic and terrestrial processes control the projected environmental conditions, we evaluate model sensitivity to a suite of climate and land-cover change (Thorn et al. 2017).

Upper Merrimack River watershed
The UMRW is located in south-central New Hampshire, USA.The UMRW comprises 8,000 km² of New Hampshire and drains through a point just south of the city of Manchester, New Hampshire (Latitude: +43.6575,Longitude: -71.5005;Fig. 1).The watershed is currently home to 410,000 people, with ongoing population growth and land-cover change through extensive residential development, and is also an important tourist destination.Annual precipitation currently averages 1,100 mm yr -1 and is evenly distributed throughout the year.Mean annual temperature is currently 8.2°C, with mean annual minimums and maximums of -8.3°C and 22.5°C, respectively.Land cover consists of a mix of deciduous and evergreen forest (82%), urban (4.2%), https://www.ecologyandsociety.org/vol22/iss4/art18/agriculture (3.8%), wetland (5.9%), and open water (4.2%).As the dominant land cover, forests are a critical influence on the aquatic environment, including water supply, water quality, and aquatic habitat.The watershed sits at the boundary between strong and weak winters, and is expected to transition to greater and more variable precipitation and warmer temperatures in coming decades (Hayhoe et al. 2007, Burakowski et al. 2008, Wake et al. 2014).The area is also experiencing rapid population growth and ongoing residential development, which is expected to continue into the future (U.S. Environmental Protection Agency (USEPA) 2010).The Merrimack watershed is ranked in the top five watersheds nationally in terms of projected changes in water quality due to increase housing density on private forest lands and ranks at the top of the list in terms of private forest land projected to experience increased housing density in the USA (Stein et al. 2009).These changes are leading to increased water use, nitrogen (N) fluxes, chloride concentrations, and other changes.Thus, this watershed is an ideal location to focus efforts in understanding ecosystem service change.

Environmental indicators
Identifying clear and tangible environmental indicators is a necessary step to communicate the response of complicated human-environmental systems to change (Müller and Burkhard 2012) and to track the performance of environmental programs, regulations, and agencies (Boyd 2004).A multidisciplinary group consisting of ecologists, hydrologists, engineers, economists, and decision scientists developed a set of potential environmental indicators important to the study region.We used an iterative and collaborative process to distill a set of climate, terrestrial, and aquatic ecosystem variables to environmental indicators.Changes in these indicators in response to different climate and land-cover scenarios were then assigned a relative value by participants in a set of multicriteria decision workshops described elsewhere (Mavrommati et al. 2017, Murphy et al. 2017).
The indicators represented multidecadal average conditions in the UMRW for the present  and projected for the end of the 21st century .Candidate indicators had relevance to climate, land, and water domains.Terrestrial indicators were based on land-cover scenarios (Thorn et al. 2017) and the Northern Research Station Climate Change Atlas (Iverson et al. 2008).Climate indicators were based the Geophysical Fluid Dynamics Laboratory (GFDL) CM2.1 global climate model (GCM) simulations using the SRES A1Fi (higher) and B1 (lower) emissions scenarios (Nakicenvoic and Swart 2000) that were statistically downscaled using the asynchronous regional regression model (Stoner et al. 2012, Wake et al. 2014).The downscaled GFDL simulations were used as they provide a reasonable representation of climate across the northeast USA (Hayhoe et al. 2007).Given the multiple land-cover and emissions scenarios used in this study, we only had the capacity to use output from one GCM.We used downscaled global climate model simulations for the Franklin, New Hampshire meteorological station, centrally located within the study domain.Water-related environmental factors were derived at the watershed scale using the coupled terrestrial-aquatic process-based model to predict ecosystem function consistent with the land-cover and climate scenarios.Generating indicators based on model simulations at daily time-steps permits the clear definition from key variables such as concentrations, temperature, or flow volumes compared with key thresholds relevant at the scale of days.Translation of environmental indicators to indicators of ecosystem services and valuation of these services is discussed later in the special issue (Mavrommati et al. 2017).

Coupling of terrestrial and aquatic models
We coupled the existing forest (PnET-CN) and aquatic (FrAMES) models to simulate hydrological and water quality characteristics related to ecosystem services, consistent with forest responses to climate, at regional scales.The detailed description of the individual PnET-CN and FrAMES models are presented in Appendix 1, sections 1.2 and 1.3.The PnET-CN model is widely used to simulate forest water, carbon (C), and N dynamics (Aber et al. 2005, Ollinger et al. 2008).The FrAMES model (Wollheim et al. 2008a, b, Wisser et al. 2010, Stewart et al. 2011, 2013) is a spatially distributed gridded river network model that has been applied extensively at various spatial scales ranging from watershed to global domain for simulating N dynamics in rivers (Wollheim et al. 2008a, b, Stewart et al. 2011), runoff/discharge dynamics (Vörösmarty et al. 1998, Wisser et al. 2010), river water temperature (Stewart et al. 2013), and chloride concentration (Zuidema, Wollheim, Mineau, et al., unpublished manuscript).
Typically, FrAMES has a land-surface hydrology component that operates independently of forest dynamics.Here, we substitute https://www.ecologyandsociety.org/vol22/iss4/art18/PnET-CN predictions of runoff and N to load material from forests to river networks.
The PnET-CN model accounts for the influence of photosynthesis on evapotranspiration and nutrient uptake, forest age, and plant physiological responses (stomatal conductance) to changing CO 2 (Ollinger et al. 2008).Together, these factors control C sequestration, nutrient export, and runoff generation immediately relevant to several ecosystem services.In the coupled model (Fig. 2), PnET-CN calculates daily runoff and dissolved inorganic nitrogen (DIN) flux from forest rooting zones.These outputs are then partitioned into shallow groundwater or surface (quick flow) flow paths, with different characteristic travel times.In urban regions, precipitation and snowmelt on hydrologically connected impervious areas run directly to the stream network, with the remainder infiltrating to lawn areas.Chloride, a potential stressor of aquatic biota, from snowmelt on the road-salt-treated fraction of impervious areas is transported conservatively following the soil and groundwater flow paths (Zuidema, Wollheim, Mineau, et al., unpublished manuscript).To link with the aquatic network, we also incorporated the role of terrestrial flow paths and riparian zones in regulating DIN loads.The PnET-CN model predicts leaching from the forest rooting zone.To account for retention along terrestrial-riparian flow paths, we applied a constant retention factor of 70% to all leachate, consistent with reactivity of riparian zones (Green et al. 2009) or buffers that average about 25 m (Mayer et al. 2005).We account for DIN loading from urban and agricultural areas using an empirical relationship between DIN, land use, and flow found in other New England watersheds (Wollheim et al. 2008a, b, Stewart et al. 2011, Mineau et al. 2015).Water temperature in terrestrial runoff is modeled as described in Stewart et al. (2013).Water temperature and DIN inputs from land are further modified by instream processes as water flow through the river network.FrAMES routes water, DIN, water temperature, and chloride originating from forests and urban and agricultural land through the river network, accounting for additional point sources of contaminants, such as N in waste water treatment plant (WWTP) effluent and thermal loads from power plant cooling water processes (Stewart et al. 2013, Miara andVörösmarty 2013).It also accounts for instream transformations during routing, such as dilution, denitrification (Wollheim et al. 2008b, Stewart et al. 2011), and instream temperature re-equilibration (Stewart et al. 2013).River discharge, solute masses, and thermal fluxes are propagated through the gridded river network using a linear reservoir routing scheme.Thus, PnET-FrAMES accounts for terrestrial-aquatic and upstream-downstream linkages to provide a mutually consistent climate, land, and water response in future scenarios.The model was first applied to historical conditions (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) to test against field observations using downscaled reanalysis meteorological observations and current land cover.Meteorological forcing for future scenarios used statistically downscaled climate data from Hayhoe et al. (2007), and land cover from Thorn et al. (2017).A detailed summary of the model input data for contemporary and future scenarios is described in Appendix 2.

Model parameterization
Model parameters were specified a priori.For future scenarios, parameters that are responsive to human management were altered to be consistent with each scenario narrative (Table 1).Key parameters include the fraction of hydrological connected imperviousness (affecting flow regime, total runoff), road salt application rates (affecting water quality), suburban and agricultural DIN nonpoint loading, and WWTP DIN point loading (affecting N exports).We are using the term "land cover" in describing future changes to the landscape.Management decisions that are consistent with the land-cover narratives (e.g., stormwater and wastewater infrastructure) are included in the parameterization of FrAMES-PnET.Specific parameterizations and how they vary by scenario are provided in Appendix 3.

Model testing with observed data
To validate the linked suite of models, we compared model simulations with data from headwater and downstream gauge locations to determine how well the spatial and temporal distribution of inputs, dilution, and processing are simulated.Observed data represented aquatic conditions at a total of 106 stations located throughout the Merrimack River and the neighboring Piscataqua River watersheds (Fig. 1).The latter is included in contemporary simulations to increase the pool of observations.Mean daily discharge data  were available from the U.S. Geological Survey (USGS) (n = 41 stations).Additional observed variables, including water temperature and concentrations of DIN and chloride came from two networks of in situ sensors: the High Intensity Aquatic Network (HIAN) (n = 12 stations) and Lotic Volunteer Temperature Electrical Conductivity and Stage (LoVoTECS) network (n = 53 stations).Data from the HIAN (Mulukutla et al. 2015) and LoVoTECS network (LoVoTECS 2012) are available through the Data Discovery Center (http://ddc.unh.edu) over the timeframe of 2012-2014 and are also reported in Contosta et al. (2017).We compared predictions to observations of both mean annual and mean summer conditions using contemporary land-cover and climate drivers.We calculated model prediction errors (MPE) and root mean square errors (RMSE) for temperature, runoff, and nitrate and chloride concentrations.We compared predictions in headwater streams (to test predictions of inputs from land) and along the Merrimack River main stem to test the regionalization and aquatic process component of the model.

Choice of environmental indicators
The iterative and collaborative process of environmental indicator selection identified ten indicators, from an initial set of 43 (Table A2.1), to represent a comprehensive suite of environmental conditions for the UMRW.The experts group used the criteria that indicators should hold perceived relevance to the general public, should equitably represent the three domains, and should be constrained by existing modeling capability.The selected indicators represented three indicators for each of land and climate domains, and four indicators for the water domain.The indicators are defined in Table 2 and explained in detail in the supplementary material (Appendix 2).The ten selected indicators are similar to several included in comprehensive lists generated previously (de Groot et al. 2010, Burkhard et al. 2012).However, our final indicators differed from previous lists to be more relevant to local residents.For instance, we used total forest cover as a percentage of the watershed instead of wood biomass stock in units of mass per area (de Groot e al. 2010) based on lack of relevance of wood biomass to the general public living in the watershed.The water indicators emphasized basin scale and subannual estimates of water conditions, requiring space-and time-varying aquatic modeling.The choice of indicators in representing specific ecosystem services is discussed more thoroughly elsewhere in this special issue (Mavrommati et al. 2017).

Regional validation of the coupled model
The component ecosystem models used in this analysis have been previously validated individually.The PnET-CN model was validated for water and nutrient balance, (Ollinger et al. 2002, 2008, Aber et al. 2005), whereas FrAMES was validated for river discharge using a simple water balance model (Vörösmarty et al. 1998, Wisser et al. 2010), water temperature (Stewart et al. 2013), chloride (Zuidema, Wollheim, Mineau, et al., unpublished manuscript), and riverine DIN (Wollheim et al. 2008a, b, Stewart et al. 2011).In this study, terrestrial vegetation processes in PnET-FrAMES drive regional runoff, nutrient, chloride, and water temperature loading that is then routed through the network.
The PnET-FrAMES model captures the range of variability in mean annual DIN concentration loaded via nonpoint inputs in headwater streams of different land cover, which reflect catchment loading prior to river processing.Predicted DIN loading concentrations are higher than observed for forested catchments (n = 4), but show little bias for developed catchments (n = 6) (Fig. 3; Figs.A4.1c, A4.2c).The PnET-CN model is known to overestimate forest DIN (Zhou, Ollinger, Glidden, et. al., unpublished manuscript), and continues to do so here, even though we added a flow path/riparian retention term to PnET-FrAMES.As a result of high bias in forested sites, the percent errors for DIN in Fig. 3 are higher than other variables.However, absolute error at these sites is very low compared with the more developed sites (Figs.A4.1c, A4.2c).Overall bias is low in summer (median MPE = -3.9%,IQR = −53.0-111.4%)but high annually (median MPE = 17.5%,IQR = −41.8-186.4%).Dissolved inorganic N concentrations in headwater streams reflect loading because aquatic net removal has not had the opportunity to impact concentrations due to short surface water flow paths in small catchments (Wollheim et al. 2006).Observations along the river main stem (basin profile, Fig. A4.3a) indicate the model reasonably represents regional loading, dilution, and aquatic N removal throughout the river network (Fig. A4.4a).

Regional framework for land, climate, and water indicators to 2100
We characterized the suite of environmental responses to two land-cover and two climate scenarios selected for the valuation workshop (Mavrommati et al. 2017) meant to represent extremes of potential change in future conditions.We drove the model with the Backyard and Small Community Food Land-Cover scenarios (Table 1), with higher (A1Fi) and lower (B1) climate scenarios, respectively (Backyard/High emissions and Small Community Food/Low emissions).
Changes in the land indicators (Table 2; Agricultural Cover per person and Forest Cover) directly reflect the land-cover scenarios (Thorn et al. 2017) and have important direct ecosystem service implications.Briefly, in the Backyard scenario, watershed forest cover declines from 80% in 2010 to 60% in 2100, to accommodate increased rural and suburban development (Fig. 4).In the Small Community Food/Low scenario, forest declines to 64% in 2100, replaced by agricultural lands, primary for pasture and hay.By 2100, population increases by 170% in the Backyard/High scenario and declines by 5% in the Small Community Food/Low scenario.The increase in agricultural land relative to population under the Small Community Food/Low scenario (~1.0 acres per person in 2100 compared with 0.2 in 2010) increases the potential for per capita local food production compared with the contemporary period (Thorn et al. 2017).The Maple Presence indicator, which primarily represents the geographic distribution of suitable habitat for Acer saccharum (Iverson et al. 2008), declines from 49% of total forest cover in the present day to 27% and 31% in the Backyard/High emission and Small Community Food/Low emission scenarios, respectively.The declines result from shifts in climate associated with each land-cover scenario.Projected changes in climate indicators in the UMRW are strongly determined by emission scenario.Present-day climate trends (Hodgkins and Dudley 2006, Hayhoe et al. 2007, Burakowski et al. 2008) continue through 2050 regardless of emission scenario (Fig. 5).Around 2050, the rate of warming is projected to increase rapidly under the Backyard/High emission scenario.With the high emission scenario, the hot days indicator increases sharply, averaging >45 d above 32°C after 2075 (Fig. 5a).Under the lower https://www.ecologyandsociety.org/vol22/iss4/art18/emission scenario, the multiyear mean for hot days is considerably lower (<17 d yr -1 ) (Fig. 5a).The Snow Days indicator declines under the higher emissions (averaging 18 d yr -1 ) compared with contemporary mean (60 d yr -1 ) and relative to the lower emissions scenario (47 d yr -1 ) (Fig. 5c).Again, the change is accelerated after 2040.The Recreation Days indicator changes little under both the low and high emission scenarios (Fig. 5b).Minimal changes in recreation days results from loss of recreation days due to days becoming too hot during the summer being offset by additional days in the shoulder seasons (spring and fall) that were previously too cool falling in the recreation day range.Environmental indicators in the water domain reflect the combined influence of projected climate and land-cover change, as well as terrestrial and aquatic processes.In general, all water indicators reflect increasing degradation of aquatic ecosystem services for the Backyard/High scenario, but remain relatively unchanged for the Small Community Food/Low scenario (Fig. 6).The Fish Habitat Loss indicator is projected to rise substantially in the Backyard/ High scenario to average about 40% of total stream and river length after 2075.Only slight increases, and fewer extreme events, occur under the Small Community Food/Low scenario.The Nitrogen Export indicator increases steadily over the 100-yr period for the Backyard/High scenario to 995 tons N yr -1 up 3100% from current conditions.Under the Community/Low scenario, however, N export changes relatively less, increasing only 70% to an average of 52 tons N yr -1 .The steady increase in the Backyard/High scenario is due to elevated N loading associated with widespread land-cover change and increased population.The pattern in Community/Low occurs because in this scenario, agricultural management activities were assumed to reflect the greater concern for the environment (optimal fertilizer application rates, improved waste water treatment plant N removal), and because it reflects a lower human population (Table 1).The Water Shortfall indicator increases in the future in Backyard/ High (averaging 6×10 6 person days [pd] yr -1 after the 2075) but remains unchanged in the Community/Low scenario.The Water Shortfall indicator increases in Backyard/High despite greater runoff that occurs due to increased precipitation (1320 compared with 1270 mm yr -1 between the higher and lower emission scenarios, respectively) and lower forest evapotranspiration (ET) due to trees responding to increased atmospheric CO 2 by increased water-use efficiency (Ollinger et al. 2002).The increase occurs because the greater population (in Backyard) reflects a greater demand that at times would be unmet.
The number of people at Flood Risk is expected to increase for both scenarios, although there is considerable interannual variability, which is expected because extreme flows only occur periodically.Total population affected after 2075 by extreme flows leading to Flood Risk is higher in the Backyard/High scenario (1.2×10 5 pd) than in Small Community Food/Low scenario (4.4×10 4 pd) because the former has a larger population to potentially be impacted, has higher flows on average, and has greater storm runoff due to increased imperviousness.
Overall, the suite of indicators suggests that relevant environmental conditions remain relatively stable until 2050 regardless of climate and land-cover scenario.Beginning in 2050, there is a potential for climate and water indicators to significantly deteriorate in the high emission scenario.The timing of the change in water indicators in Backyard/High suggests primarily a climate control, as climate changes also show rapid increases at this time (Fig. 5).https://www.ecologyandsociety.org/vol22/iss4/art18/

Using coupled models and multiple scenarios to partition influence of land cover and climate on water indicators
Climate drivers govern the hydrological and biogeochemical response of catchments, but these changes are mediated by land use and land cover through impacts of impervious surfaces, plant growth, evapotranspiration, and water residence time (Ollinger et al. 2002, Raymond and Saiers 2010, Pellerin et al. 2011, Goodridge and Melack 2012).The two scenarios described above convolve both signals (i.e., both climate and land cover differ in the two scenarios).The relative importance of climate, land cover, and interactions on the water domain indicators can be determined by analyzing results from all combinations of the suite of land-cover change scenarios (Table 1) and both C emission (climate) scenarios (Fig. 7).
The relative importance of climate and land cover varies among the indicators.For Fish Habitat Loss, climate dominates the response as the high emission scenarios have greater habitat loss than the corresponding low emission scenario across all landcover scenarios (Fig. 7a).Specific land-cover scenarios modify only moderately the degree of fish habitat loss.Fish habitat loss is primarily driven by temperature impairment in headwater streams (Fig. A4.4b), resulting in a 4-to 10-times increase in river kilometers impacted in the high compared with low emission scenarios (Fig. 5a, 7a).The Backyard scenario exacerbates fish habitat changes caused by climate, whereas the three other landcover scenarios remain similar.The Fish Habitat indicator is composed of flow, water temperature, and chloride impairment (Table 2).The greater degradation in Backyard occurs because of https://www.ecologyandsociety.org/vol22/iss4/art18/ the greater abundance of impervious surfaces, with increased runoff from the land leading to greater spatial and temporal probability of crossing chloride and water temperature thresholds.Residential development and associated impervious surfaces in the Backyard scenario result in greater summertime heat flux to streams through stormwater runoff relative to other build-out scenarios.In the Community (both Small Community Food and Large Community Wild) scenarios, reduced effective imperviousness from coupling low-impact design buildout with concentrated development results in negligible increases in heat flux impacts to rivers and no additional fish habitat loss relative to current land cover and projected climates.
For Nitrogen Export, climate is also the primary determinant of indicator response, but interactions with land cover are also evident (Fig. 7b).The Nitrogen Export indicator is higher in all higher emission scenarios relative to lower emission scenarios.
Higher emissions scenarios result in greater coastal N flux than lower emissions scenarios due to greater loading to the river system.First, the greater precipitation in high relative to low emission scenario results in higher atmospheric N inputs as well as higher water runoff.Second, greater forest N leaching results from higher temperatures under future scenarios due to declining uptake by conifers.Third, the river system has reduced capacity to retain N under higher flow regimes (Wollheim et al. 2008b).Small Community Food shows little response under the low emission climate scenario (as discussed above) similar to if land cover remained unchanged.Both the Large Community Wildlands and Backyard scenarios show additional N export increases relative to other land-cover scenarios, indicating an interaction between land cover and climate.These scenarios have similar high populations compared with Small Community Food, and result in N loading into rivers via domestic waste, either through septic systems (Backyard) or WWTPs (Large Community), which lead to increased waste N inputs.
The Large Community Wildlands scenario projects greater N export than the Backyard scenario (Fig. 7b), despite the former conserving more forest, and upgrading WWTPs to higher treatment levels with lower per capita human waste N inputs.The logistic loading function that is used to parameterize nonpoint N inputs to rivers as a function of land cover (Appendix 3) assumes relatively low N increases for low to moderate density urban and agricultural development (average increases to about ~21% urban + Ag land cover for Backyard and 9% for Large Community Wildlands).The mechanistic explanation for this pattern is that N removal by terrestrial ecosystems remains high up to certain thresholds of natural ecosystem loss, consistent with a number of previous studies (Groffman et al. 2004, Wollheim et. 2005).The Backyard scenario assumes a low density of future land-cover change that is often below the threshold.Thus, human waste associated with population growth, which is managed through distributed septic systems and not WWTP, is mostly retained, and there is relatively little change in N inputs.In contrast, the Large Community Wildlands scenario assumes domestic waste N is transferred to WWTP with Total Nitrogen (TN) removal efficiencies of 90%.The assumption of low response to land-cover change embedded in the loading function is derived from one location with a certain social and biogeophysical context (e.g., rural/suburban septic infrastructure, wetland abundance, etc. Wollheim et al. 2008b).Although the generality of this loading function to the entire domain should be tested more thoroughly, the validation results suggest it is reasonable and therefore useful for exploring potential future responses.
An additional small driver of the difference in N export between Backyard and Large Community Wildlands was that the former also presents greater opportunity for instream removal than the latter because loading occurs to smaller streams (as nonpoint sources) compared with waste water effluent to larger rivers (as point sources).The resulting longer lotic travel time and greater exposure to benthic surfaces provide more opportunities for N removal through denitrification (Mineau et al. 2015).However, in our simulations, only an additional 1.6% of total network inorganic N inputs was removed (denitrified) in the aquatic system under the Backyard scenario compared with Large Community Wildlands.
Water shortfalls are controlled more by land-cover and population scenarios than by climate scenarios (Fig. 7c).The lack of a strong climate driver occurs because conditions are expected to become increasingly wet under both lower and higher emissions scenarios, alleviating water shortages.The highly dispersed development of the Backyard scenario places an increasing population throughout headwater catchments, dependent on well water.In contrast, the Community scenario places more people on public water supplies that obtain their water from larger catchments.The Water Shortfall indicator we used relates water demand accumulated through the river network (a function of population and per capita demand) with surface water discharge in the grid cell.In the Backyard scenario, people are presumed to depend primarily on wells, not urban water supply infrastructure, where small catchments have limited hydrological storage.Thus, demand more frequently outstrips supply in these local areas.
Runoff increases in the higher emission scenario compared with the lower emission scenario, which reflects a combination of greater annual precipitation predicted under a warmer climate and greater forest water-use efficiency under a C-enriched atmosphere (Fig. A4.4c, d).As a result, the high emission scenario actually reduces water shortfalls relative to the low emission scenario (Fig. 7c).
Flood risk is projected to increase under both emission scenarios (Fig. 7d), resulting from greater precipitation (Fig. A4.4c) occurring during extreme events (Hayhoe et al. 2007, Trenberth 2011).Backyard and Large Community Wildlands scenarios project a larger flood-affected population due simply to the larger populations compared with other scenarios.Generally, the higher emission scenario has greater flood risk than the lower emission scenario, with the exception of the Backyard scenario.High emissions are expected to have slightly higher precipitation and runoff relative to low emissions.The Backyard scenario differs in that more people live further north in the watershed than in the Community scenario (concentrated on existing urban areas, Fig. 1), which may expose these populations to more extreme flows.

Spatial variability of fish habitat loss
The spatially distributed modeling approach used here is able to identify changes in environmental indicators relative to human settlement.As an example, the Fish Habitat Loss indicator, as the fraction of years with unsuitable habitat, is greatest in the southern part of the study domain (Fig. 8).Greater impairment in the southern portion of the watershed occurs because warmer air temperatures increase water temperatures as rivers flow from north to south (Stewart et al. 2011) and because denser road networks relative to the northern portions of the watershed experience greater winter road-salt application.
The Contemporary period shows elevated fish habitat loss in the smaller streams in the southern half of the watershed, but not in the larger rivers (Fig. 8a).Meteorological droughts in 1981 and 1993 are responsible for the contemporary fish habitat loss, which was caused by the low flow component of the indicator.It is possible that future climate from the GFDL GCM used in these scenarios does not capture longer periods of summer dry conditions later in the 20th century (Hayhoe et al. 2007, Wake et al. 2014).Thus, the scenarios presented here may underestimate the effect of the low-flow criteria on future fish habitat loss (Fig. 8).
Future scenarios are likely to intensify fish habitat loss near population centers in the high emission scenarios (Fig. 8b, d, top), while spreading through more of the previously unimpaired part of the watershed.In contrast, the low emission scenarios tend to shrink the footprint of the Fish Habitat Loss indicator, although retaining the likelihood of impairment near the population centers (Fig. 8c, e, bottom).The Backyard scenario shows a stronger response than the Community scenario because greater impervious cover leads to both increased road-salt loading and increased temperatures from increased summertime impervious surface runoff.
The greater Fish Habitat Loss indicator in large rivers compared with the contemporary period is particularly evident across all future scenarios, resulting in large part because of warmer water temperatures.Large river water temperature responds disproportionately because atmospheric warming further warms water as it flows downstream during summers, whereas headwater streams remain cooler because of groundwater recharged during cooler times of year.The spatially distributed, time-varying modeling approach used here allows consideration of heterogenous impacts to environmental function that can be used in valuation studies.

Role of land-cover management
Despite the predominant influence of climate on the environmental indicators in the UMRW, land-cover management can mitigate degradation to certain attributes of the ecosystem (e.g., Fig. 8).Build-out scenarios directly influence land-based environmental indicators such as Forest and Agricultural Cover.Scenarios that maintain greater forest land have slightly greater C sequestration (Fig. A4.4e), whereas scenarios with greater agriculture increase food supply (Thorn et al. 2017) and resiliency (Tilman et al. 2002).In addition, these land covers affect numerous sociocultural resources, including water-related ecosystem services.Depending on the relative value placed on different indicators (Murphy et al. 2017), the landscape can be managed to maximize ecosystem services that are most important.https://www.ecologyandsociety.org/vol22/iss4/art18/Coupled terrestrial-aquatic models express the understanding of processes relevant to specific land-cover decisions, permit exploration of relevant tradeoffs (Antle andCapalbo 2001, Bennett et al. 2009), and are essential in understanding the response of ecosystem dynamics to changes in both climate and land cover.Results from PnET-FrAMES provide a measure of the potential impact of each land-cover scenario and each climate scenario on different ecosystem service indicators.As regional managers have little influence over climate change and because we are currently on a higher emission trajectory (Pachauri et al. 2014), managers should consider potential land cover in conjunction with the higher emissions climate scenario to guide future planning.Projections from PnET-FrAMES suggest that infilled development (Community) scenarios exhibit two types of paradigms with respect to specific aquatic environmental indicators under a high emission future (Fig. 7).The three landcover scenarios with limited expansion of the residential footprint (Constant, Small Community Food, Large Community Wild) show smaller levels of Fish Habitat Loss and Water Shortfalls than the dispersed buildout (Backyard).The Large Community Wildlands scenario represents a doubling of population similar to Backyard scenario, but has much lower environmental detriment based on these two indicators.However, the Large Community Wildlands scenario does show increased impact on the Nitrogen Export indicator (more point sources) and Flood Risk indicator (more people along large river corridors) relative to Backyard, suggesting that there are tradeoffs that must be further managed under an infilled paradigm.
The scenarios tested are simple and do not encompass numerous potential approaches that managers may use to mitigate the adverse effects of the investigated environmental stressors.The scenarios do not account for potential nonpoint N runoff controls including managed riparian buffers (Mayer et al. 2005) or stormwater control (Collins et al. 2010) that would be consistent with Community scenarios.Nor do the scenarios account for novel flood control measures such as blue-green infrastructure (Thorne et al. 2015) or explicit development outside of flood zones, because flood zones are mapped at resolutions finer than our simulations.Therefore, the flood risk identified here may potentially be mitigated through engineered interventions.
Finally, potential management interventions such as constructed stream refugia (Sedell et al. 1990, Isaak et al. 2015) are not considered but could mediate warming stream temperatures for instream biota, a threat under all build-out scenarios.The scenarios presented here identify those environmental stressors that are most likely to need more directed research and management to mitigate changing populations and expected warming.

CONCLUSION
Process-based terrestrial-aquatic models are essential for quantifying potential future environmental impacts due to spatially distributed changes in land cover and climate.Here, we simulated dynamic interactions across land and water domains at a daily temporal resolution and a spatial resolution that accounts for the heterogeneity and associated processes of the UMRW in New Hampshire.Our approach linked models of intermediate complexity to connect terrestrial and aquatic domains and allowed us to explore the impacts of changes in climate and land cover on aquatic ecosystems through the use of different scenario combinations.Process-based models enable understanding of macroscale ecosystem responses based on underlying ecosystem processes.
Water is arguably the most fundamental resource for society, and it is therefore critical to quantify how water quantity and quality will respond to projected changes in land cover and climate.Our results suggest that climate is expected to be the most influential driver of water-related impairment in the UMRW, although impacts to climate and aquatic environment conditions are likely to be modest until mid-century.(Wythers et al. 2013).e.g., at the temperature of 35 o C, the RA algorithms estimate a 37% reduction in foliar respiration relative to that using previous version.Wythers et al. (2013) reported that averaged across four boreal ecotone sites and three forest types at year 2100, the enhancement of NPP in response to the combination of rising CO2 and warming was 9% greater when RA algorithms were used, relative to responses using fixed respiration parameters.
Biogeochemical monitoring for 50 years at the Hubbard Brook Experimental Forest in New Hampshire has revealed N export in stream water has steadily declined and is presently just a small fraction of atmospheric N input, despite negligible changes in aboveground biomass.It implies that the forested ecosystem has shifted to a net N sink (Yanai et al. 2013), which could not explained by the previous theory.The "missing" deposited N were thought to accumulate in the mineral soil, or be lost in gaseous form.
Processes of N gas losses (i.e., N2O, NO, and N2) through nitrification and denitrification were added in PnET-CN to enhance its N cycling and to investigate the role of denitrification in the missing N sink (Zhou et al. In Prep).We used first order kinetics to estimate N gas losses and partitioned N2O, NO, and N2 based on soil water content.In this study, the parameter, denitrification constant was set to the averaged value of 0.03 (McCray et al. 2005) to represent a more general pattern for a large region, which could be potentially underestimate in mountainous areas.

FrAMES
FrAMES, the Framework for Aquatic Modeling in the Earth System, is a spatially distributed gridded river network model that has been applied extensively at various spatial scales (Wollheim et al. 2008a, 2008band Stewart et al. 2011, Vörösmarty et al. 1998, Wisser et al. 2010, Stewart et al. 2013 and Zuidema et al. In Prep).FrAMES incorporates a number of dynamically linked modules that operate on a daily time step.

Future land cover scenarios
To demonstrate the coupled model, and to develop indicators required for the ecosystem services valuation (Mavrommati et al. 2017), we focus on two land cover scenarios expected to show the largest range in changing ES, the "Backyard Amenities" (Backyard) and the "Small Community with Promotion of Local Food" (Small Community Food) (Thorn et al. this issue).Table 1 presents key differences between these land cover scenarios.The Backyard land cover scenario, which prioritizes large building lots and incurs increased transportation related energy consumption, was paired with higher greenhouse gas emission (A1Fi) scenario.Conversely, the Small Community Food scenario reduces transportation-related energy consumption and is more consistent with the lower greenhouse gas emission (B1) scenario.For subsequent analyses into the specific roles of climate and land cover change, we consider the responses of the suite of land cover (Thorn et al. this issue) and emission scenarios.
We simulate future drivers of terrestrial and aquatic ecosystems using projected impervious cover, population density, and land cover from Thorn et al. (this issue) for each future decade from 2020 to 2100.Land cover and population data was aggregated using the same methodology as for contemporary land cover.We assumed the proportion of each forest type remains constant in the future.) is then calculated as the maximum of local population or the population that is expected to have insufficient supply, except in densely populated cities (  > 350   −2 ) where supply is assumed available elsewhere in the watershed.This density is consistent with areas of known public water supply systems in the UMRW.The demand of these areas is transmitted downstream, so the high demand can trigger water shortfalls in downstream cells.

Water domain
As formulated,   considers both rural and urban water consumers in a simple consistent framework that is flexible across development patterns.The indicator assumes that streamflow is an appropriate indicator of water supply even where actual abstractions may come from groundwater resources.
Flood attenuation (  [ ]).Flood attenuation was selected as an indicator of potential impact to human and infrastructure safety.Flood attenuation represents the spatial and annual sums of population potentially affected (  ) by local flooding, times the number of days per year persons were potentially affected.To calculate affected population, we assume that the entire population of a grid-cell (at ~1.5 km2) is potentially affected when daily mean discharge in the grid cell from PnET-FrAMES exceeds an estimate of the 100-year flood discharge.The threshold defining flood runoff was spatially distributed dependent on mean April precipitation, channel slope, upstream wetland fraction, and catchment area following Olson et al. (2009).
Fish habitat ( ℎ [𝑘𝑚]).Fish habitat was selected as an indicator of potential non-use or aesthetic value, as well as recreational value of sport fishing.Fish habitat is calculated as the fraction of the total length of all streams and rivers in the watershed that exceed at least one threshold associated with either water temperature, salinity (as chloride), or instream flow at least one time in each year.We use a 7-day rolling mean water temperature and determine when this value exceeds the median of freshwater fish tolerances (29.2°C) from Eaton and Scheller (1996).We use 4-day rolling mean of chloride to determine when it exceeds the USEPA chronic water quality criteria for chloride (230  Cl  −1 ) (USEPA 1988).We use the 7-day rolling mean of discharge to determine when it is below an estimate of in-stream flow requirement equal to present day 7Q10, or the 10 th -percentile of annual minimum 7-day rolling mean discharges.We estimate present-day 7Q10 in runoff as 0.122 [  −1 ] using USGS daily discharge data for 12 stations in the Merrimack River Watershed with at least 17 years of data.
Nitrogen export (  [ N  −1 ]).Nitrogen export was selected as an indicator of the watershed regulation service of processing excess anthropogenic dissolved inorganic nitrogen (DIN) prior to export to sensitive estuarine habitats.

Fig. 1 .
Fig. 1.Merrimack and Piscataqua River networks and watershed boundaries showing distribution of land cover and validation data stations.

Fig. 2 .
Fig. 2. Conceptual diagram of the terrestrial and aquatic ecosystem model coupling and factors influencing environmental indicators.

Fig. 3 .
Fig. 3. Box plots of model percent error (MPE) across all sites (MPE = (P-O) / O * 100) for aquatic variables during annual (a) and summer (b) time periods (RO: Runoff; DIN: Dissolved Inorganic Nitrogen, Temp: Water Temperature; where, P is the predicted values by the model and O is the observed values at different sampling sites).

Fig. 4 .
Fig. 4. Time series of land environmental indicators for Backyard Amenities / High Emissions and Small Community Amenities / Low Emission scenarios through 2100 for (a) Agricultural Cover (acres/person), (b) Forest Cover (proportion), and (c) Maple Suitability (% forest cover).

Fig. 5 .
Fig. 5. Time series of climate environmental indicators for Backyard Amenities / High Emissions and Small Community Amenities / Low Emission scenarios through 2100 for (a) Very Hot days (days per year), (b) Comfortable days (days per year), and (c) Snow Cover days (days per year).

Fig. 6 .
Fig. 6.Time-series of water environmental indicators for Backyard Amenities / High Emissions and Small Community Amenities / Low Emission scenarios through 2100 for (a) Fish Habitat Loss (% river length impairment per year), (b) Nitrogen Export (10 6 x kg N per year), (c) Water Shortfalls (Million (M) person days per year), and (d) Flood Risk (Thousand (k) person days per year).

Fig. 7 .
Fig. 7. Distribution of water indicators between 2075-2100 for eight scenarios (four land use and two climate; Table1) compared with contemporary range, including (a) Fish Habitat Loss (km d yr -1 ), (b) Nitrogen Export (kg N yr -1 ), (c) Water Shortfalls (P d yr -1 ), and (d) total Flood Risk population-duration (P d).Boxes in panels (a)-(c) represent first to third quartiles of data with line defining median, whiskers extend to 1.5 x IQR, outliers as circles beyond.Bars in (d) represent the cumulative total population potentially flood affected during the period.Highlighted columns depict scenarios included in workshop evaluation(Mavrommati et al. 2017) and focused on Figs.4-6.

Fig. 8 .
Fig. 8. Map showing the spatial variability of one of the critical ecosystem service indicators (Fish Habitat Loss) as proportion of years when fish habitat loss exceeds impaired levels for the study domain under (a) Contemporary (1980-2005), (b) Backyard Amenities (2075-2100) high emission (top) and low emission (bottom), and (c) Community Amenities (2075-2100) high emission (top) and low emission (bottom).Location of urban population (density > 200/km²) is shown in red in each scenario.

Four
indicators were chosen for the water domain that represent a diversity of ecosystem services.The four indicators were derived from coupled PnET-FrAMES model output forced using the same land cover and climate datasets used to derive indicators for atmosphere and land domains.The four water domain indicators were expressed in terms of potential stressors to relevant ecosystem services.Associated supply and demand of ecosystem services related to water also evolve with urbanization(Wollheim et al. 2015).Water supply (  [ ]).Water supply was chosen as an indicator of potential stress to human provisioning needs.Water supply represents the spatial and annual sums of population potentially affected (  ) by limited water availability, times the number of days per year persons were potentially affected.To calculate potentially affected population, we first use to calculate the daily flow accumulation (FA)(Tarboton et al. 2009) at each grid cell to estimate water availability for the population in the grid cell ().We then accumulate the downstream difference between water supply ( −   ) (stream flow minus environmental flow) and total human consumptive demand ( [  −1 ] =  [] * [  −1  −1 ]) where  is a per capita water demand of residents in New Hampshire, which 284 L/p/d(Horn et al. 2008). = (max( −   , 0) − )

Table 1 .
Thorn et al. 2017ls of different land-cover scenarios used in this study (based onThorn et al. 2017).

Table 2 .
Definitions and units of selected environmental indicators.(P = Person, d = day, yr = year)

Table 2
Modeled water temperature across gauging stations is low for mean annual temperature (median MPE = -17.5%,IQR=−32.7-−9.2%)withlittlebias during summer (median MPE = -2.3%,IQR=−8.1-5.5%)(Fig.3b).Because the model does well during summer when temperatures are high, the water temperature component of the Aquatic Habitat indicator, which is mostly driven by summer conditions, is robust.No differences in bias were seen across different river sizes.
After 2050, aquatic impacts are projected to increase rapidly in the high emission scenario, and degree of change is further influenced by land cover.Because recent analysis points to the likelihood of a high emissions scenario(Pachauri et al. 2014), the region should prepare for changes in ecosystem services.Land-cover management is the leverage point that local and state decision makers have to adapt to global climate change, and management interventions that mitigate changes in population and temperature will need to move beyond simple land-cover/land-use formulations analyzed here.
(Stoner et al. 2013 Anthropogenic Chloride Loading (NACL) module(Zuidema et al.In Prep).Typically, FrAMES has a land surface hydrology component that operates independently of forest dynamics.In this study, FrAMES has been coupled with PnET to simulate the dynamic interactions among terrestrial and aquatic processes in the UMRW across a spectrum of land cover and climate projections for the region.Here, we substitute PnET-CN predictions of runoff and nitrogen to load material from forests to river networks.In the PnET-FrAMES coupling, PnET-CN predicts runoff and DIN leaching from forests, whereas FrAMES simulates inputs of specific conductivity, thermal loads, and DIN loads from urban and agricultural areas.Water from the soil root zone in PnET as runoff is partitioned to surface and groundwater runoff generating pools in FrAMES, which introduce a lag in delivery to the stream network.PnET nitrogen leachate is applied to the daily runoff volume of the linked model.PnET does not consider riparian nitrogen removal, which is parameterized in FrAMES as a zero-order removal process that eliminates 75% of the forest leachate prior to entering the stream network.This value was determined by calibration to extensive headwater concentration data (Appendix 3) and deserves additional investigation.FrAMES propagates discharge and all solute loads downstream using a cascade routing method at a daily time step.statisticaldownscalingtechnique(Stoner et al. 2013, Wake et al. 2014) these data were unavailable as gridded data for our model domain.However, a check of projections at the specific location indicate they are consistent.

Table A2
Forest type (  ).In addition to services provided by forest cover, the type of forest, particularly the proportional presence of maple/beech/birch (MBB) is considered Farm land (  ).The farm land indicator was calculated directly from build-out scenario data including modeled agricultural land cover and total population within the UMRW.Farm land was expressed as the number of acres of watershed agricultural-land (  ) per resident of the UMRW ().