Impact of climate change and climate anomalies on hydrologic and biogeochemical processes in an agricultural catchment of the Chesapeake Bay watershed, USA

• Climate change resulting in increased nutrient cycling, particularly for nitro-


Introduction
Climate change has the potential to impact hydrology, nutrient cycling, and greenhouse gas (GHG) emissions from agroecosystems in the Mid-Atlantic region of United States (Huntington, 2003;Johnson et al., 2012;Najjar et al., 2009;Najjar et al., 2010;Neff et al., 2000). Climate predictions for the Mid-Atlantic suggest that precipitation (especially during winter and spring) quantity and intensity will increase, which could increase nitrogen (N), phosphorus (P) and sediment export from watersheds (Chang et al., 2001;Cousino et al., 2015). Changes in soil moisture, temperature, and nutrient cycling rates also influence GHG emissions (Maag and Vinther, 1996). Excess nutrient export from the landscape to aquatic ecosystems accelerates eutrophication and harmful algal blooms (Burgin and Hamilton, 2007), leading to undesirable changes in ecosystem structure and function (Smith et al., 1999), such as coastal dead zones (areas of low oxygen concentrations) (Diaz and Rosenberg, 2008). The Chesapeake Bay is a eutrophic estuary that has been the subject of precedent-setting regulations to limit point and non-point source pollution from its 166,000 km 2 watershed. Commitments have been made by all six Chesapeake Bay states and the District of Colombia to reduce contributions to Bay nutrient and sediment loading over the next several decades, however, a major uncertainty with Chesapeake Bay watershed mitigation strategies is the complicating influence of climate change and climate variability. There is also uncertainty over the effects of climate change on GHG emissions, which are often tied to the same hydrologic and nutrient cycles that are vulnerable to climatic shifts.
Nitrogen, P, and sediment export from the landscape to surface waters and GHG emissions are controlled by a combination of key biogeochemical and hydrologic processes. Changes in precipitation and temperature alter the timing and magnitude of runoff, soil moisture, and biogeochemical cycles (Gleick, 1989). For instance, N mineralization, nitrification and denitrification, are to a large extent controlled by factors that climate change influences, such as soil temperature and soil moisture (Butterbach-Bahl and Dannenmann, 2011). Similarly, increased soil temperatures and moisture content can influence the sorption and desorption of P, as well as immobilization and mineralization rates, all factors affecting P export (Sheppard and Racz, 1984). It is also well established that sediment transport is affected by soil moisture (Wiggs et al., 2004) and by precipitation amount and intensity (Römkens et al., 2002).
A fundamental understanding of these coupled processes under a changing climate is critical to managing N, P, and sediment export from ecosystems to sensitive coastal zones. Development of effective landscape management strategies to improve water quality requires an understanding of how processes that regulate nutrient/sediment production on the landscape are coupled with hydrologic transport to water bodies. Of particular interest is the impact of climate change on hydrologically active areas of the landscape that contribute disproportionally to watershed nutrient export (e.g., Critical Source Areas, CSAs), where active hydrologic transport and high nutrient availability coincide (Groffman et al., 2009). A better understanding of climate change and its potential influence on landscape biogeochemistry can be used to develop new strategies for protecting coastal waters and their contributing watersheds from pollution. For instance, it is entirely possible that climate change would enhance some natural ecosystem services that protect water quality. One example is through a potential change in denitrification, a natural process that transforms dissolved nitrate (NO 3 -) into nitrogen gases, and returns it to the atmosphere. Hydrologically active areas in the landscape prone to soil saturation are recognized as denitrification hotspots (McClain et al., 2003;Vidon et al., 2010) and understanding where these areas are, or will be, can help develop management approaches that utilize ecosystem services to protect water quality.
This study assesses the effects of climate change and weather anomalies on hydrology, nutrient and sediment export, and N 2 O and N 2 emissions in an agricultural sub-catchment of the Chesapeake Bay watershed located in east-central Pennsylvania. Long-term watershed monitoring data are used along with outputs from seven downscaled and bias -corrected regional climate models (RCMs) to understand the implications of climate change on hydrology, water quality, nutrient cycling, and GHG emissions, using the Soil and Water Assessment Tool-Variable Source Area (SWAT-VSA) model (Easton et al., 2008) modified to include GHG emissions (Wagena et al., 2017), and improved P cycling . We analyze the results of each climate model, as well as the ensemble model mean, for the impact on hydrology, water quality, N 2 O and N 2 emissions and nutrient cycling. Results highlight processes, places, and times that are vulnerable to climate impacts and point to watershed mitigation and adaptation strategies to address the long-term impacts of climate change.

Watershed description
The WE-38 Experimental Watershed is a sub-watershed of Mahantango Creek in east-central Pennsylvania, which drains to the Susquehanna River (Fig. 1). The watershed has an area of 7.3 km 2 and has been extensively studied as a USDA-ARS experimental watershed beginning in 1966 and monitored through the Sustaining the Earth's Watersheds Agricultural Research Data System (STEWARDS) project . The climate of WE-38 is temperate humid, with a mean temperature of 10.1°C, mean precipitation of 1080 mm yr −1 , and mean streamflow equal to 46% of precipitation Lu et al., 2015). Elevation ranges from 220 to 510 m and the land use of the watershed consists of agriculture (44.5%), forest (33.8%), and pasture (3.5%) . As part of the Appalachian Valley and Ridge Province, the watershed is underlain by fine sand and siltstones, supporting soils prone to variable source area (VSA) hydrology. Specifically, uplands feature well-drained soils with high infiltration rates while the lower landscape positions are more poorly drained soil with features that seasonally perch water and result in runoff generation by saturation excess processes (Lu et al., 2015).

SWAT model description
The SWAT model is a physically based, semi-distributed watershed model developed to assess the impact of land management practices on water availability and water quality (Arnold et al., 1998). SWAT requires meteorological (precipitation, min and max temperature, solar radiation, wind speed, and humidity), soil, land cover, and land management data to simulate surface and subsurface hydrology and various chemical, nutrient, and sediment fluxes. SWAT-VSA re-conceptualizes SWAT to account for areas of the landscape subject to variable saturation dynamics (Easton et al., 2008). In SWAT-VSA the area of each hydrological Response Unit (HRU) is defined by the coincidence of land use and wetness index class determined from a Topographic Index (TI) to differentiate areas of the landscape with respect to their moisture storage and saturation/runoff potential (Easton et al., 2008). SWAT-VSA has been shown to provide better predictions of soil moisture, runoff generation, and nutrient export than the standard SWAT model in WE-38 and similar watersheds (Easton et al., 2008). Wagena et al. (2017) further modified SWAT-VSA to predict N 2 O and N 2 emissions from agroecosystems based on soil physical and chemical processes and environmental drivers. The model uses reduction functions to model the total denitrification (N 2 + N 2 O) and nitrification rates, and then applies a ratio method to partition N 2 from N 2 O. Greater detail can be found in Wagena et al. (2017). Collick et al. (2016) modified the P cycling routines in SWAT to better predict the fate and cycling of surface applied P (fertilizers and manures).

Watershed model initialization and input data
SWAT-VSA was initialized with a 10-m resolution digital elevation model resampled from 0.5 m LiDAR data obtained from Canaan Valley Institute (2007) using ArcSWAT 2012 and TopoSWAT (available from https://doi.org/10.6084/m9.figshare.1342823) developed by Fuka et al., 2016. TopoSWAT automates the SWAT-VSA initialization process by creating the TI data and then overlaying the soil and TI data to develop the required database for model initialization. Soils data utilized by Topo SWAT are based on the Food and Agriculture Organization (FAO) soils database (FAO, 2007). This methodology downscales the FAO soils data, distributes the soil properties across TI classes, and has been shown to provide more accurate representation of soil properties than SSURGO . The land use characterization of WE-38 were derived from previous studies (Buda et al., 2013;Buda et al., 2009;Gburek et al., 2002;Gburek et al., 2006;Needelman et al., 2004;Veith et al., 2008) and used to parameterize the model at the field scale. There are two options to incorporate agricultural management practices in SWAT; they can be triggered by thresholds in the model (e.g. growing degree day, heat units), or they can be input into the management file as individual practices, should sufficient data exist. For this study, the later was utilized.
Multi-year agricultural field management schedules of individual fields were built into the model by Collick et al. (2015). The management information has been consistently collected for each farm in WE-38, and includes annual farmer interviews, informal updates throughout the year, and regular field observations. In order to incorporate this management information into SWAT-VSA, georeferenced field boundaries are obtained annually by GPS, reconciled against publically available aerial satellite photos, and combined to represent the changing shape and size of the fields over the 12 years of the study . The model was initialized by weather including precipitation, temperature (min and max), relative humidity, wind speed and solar radiation from 1987 to 2010 from land based stations in WE-38.

Model calibration and evaluation
The WE-38 SWAT-VSA model was calibrated and then evaluated using SWAT-CUP (SWAT Calibration and Uncertainty Procedure) (Arnold et al., 2012) using the SUFI2 (Sequential Uncertainty Fitting) optimization algorithms with the objective function set to the Nash Sutcliffe Efficiency coefficient (NSE). The SWAT-VSA model performance was evaluated based on three metrics, percent bias (PBIAS), coefficient of determination (R 2 ), and the NSE, against the historical measured data from 1989 to 1999 for model calibration and 1999 to 2007 for model evaluation. The NSE (Eq. (1)) is an indicator of the predictive power of the model and ranges from -∞ to 1 where 1 is a perfect fit between model and observed, and an NSE of 0 indicates the observed mean provides a better estimate of the data than the model (Krause et al., 2005).
The PBIAS is a statistical metric providing an estimate of over prediction (PBIAS N 0) or under prediction (PBIAS b 0) of the model based on Eq. (2).
where P i is the predicted value and O i is the observed value at time i. O is mean of the individual observations of O i , and n is the number of observations. The model was calibrated for flow, mineral phosphorus (PO 4 -), and NO 3 -using split sample approach from 1989 to 1998 and model predictive power was evaluated from 1999 to 2007 on a daily basis against measurements made at the WE-38 watershed outlet. Measured flow data consist of 15-min frequency measurements converted to a time weighted daily average. Nitrate and DP measurements are made less frequently and include both baseflow and stormflow samples, converted to a time-weighted daily average    Easton, 2018). Since the observed sediment data were only collected under baseflow conditions, and thus fail to capture storm events, calibration was not possible, we do however report model fit during the evaluation period. The predicted N 2 O emissions were evaluated (the parameters affecting N 2 O were not calibrated) against measurements made in the watershed from May 2012 to Oct 2013 (Benjamin Rau, personnel comm). The frequency of measurements varied from every four days to less than once per month depending on precipitation and fertilization events. Soil N 2 O flux was measured using vented static chambers (Parkin et al., 2003). Gas samples were analyzed for N 2 O concentration using a Varian CP3800 gas chromatograph equipped with a 63 Ni electron capture detector (Varian, Walnut Creek, CA, USA).

Climate scenario incorporation
Results from seven climate models were obtained from the North American Climate Change Assessment Program (NARCCAP) model dataset (Mearns et al., 2009). The models use dynamical downscaling -nesting a regional climate model (RCM) within a global climate model (GCM)-and provide data at high temporal (3-hourly) and spatial (50-km) resolutions, better capturing local processes (Rummukainen, 2010). Greenhouse gas concentrations in the NARCCAP future simulations are obtained from the medium-high SRES A2 emissions scenario (Nakicenovic et al., 2000). The models used include the CRCM-CCSM (Canadian Regional Climate Model, Community Climate System Model), CRCM-CGCM3 (Canadian Regional Cli- Dynamical downscaling provides some advantages over the main alternative, statistical downscaling, which uses transfer functions based on historical observations in order to relate climate model output to the climate variables of interest (e.g., Abatzoglou and Brown, 2012). Dynamical downscaling is based on the same general physical principles that underlie GCMs, such as the conservation of mass, momentum, and energy. Thus, dynamical downscaling is not tied to the historical climate parameter space that statistical downscaling approaches are tuned to. This is of particular importance for large projected changes in climate, which may be out of the range of the historical record. However, statistical downscaling has the advantage of higher resolution (due to being less computationally expensive than dynamical downscaling) and, in general, a lack of bias in simulations of the historical record as a result of the aforementioned tuning. We solve the bias problem in the NARCCAP product by using a cumulative distribution function matching method (Dixon et al., 2016;Li et al., 2010;Teutschbein and Seibert, 2013;Wang and Chen, 2014) with data from Phase 2 of the North American Land Data Assimilation System (NLDAS-2) (Xia et al., 2012). In total, three separate data sets are needed to complete the procedure: historical data as a downscaling reference, historical data from the climate model source, and future (i.e., predicted) data from the climate model source. The matching functions defined by each quantile of historical observed reference data and historical modeled climate data are applied to the future modeled climate data. The accuracy of the historical, bias corrected NARCCAP data from each of the models was verified by employing the Equiratio Cumulative Distribution Function (ECDF) matching method (Li et al., 2010;Wang and Chen, 2014) against the historical climate dataset. After verifying that the bias corrected climate models statistically reproduced the observed historical climate, the calibrated SWAT-VSA model then was forced by output of the RCMs, including precipitation, minimum and maximum temperature, solar radiation, relative humidity, and wind speed for historical model runs  and future scenarios . Note that given the small size of the watershed (7.3 km 2 ) only one NARCCAP pixel is used in the analysis.
One shortcoming of the NARRCAP product is that its RCMs are not forced by the most recent set of GCMs. NARCCAP uses GCM output from Phase 3 of the Coupled Model Intercomparison Project (CMIP3), whereas recent statistical downscaling products (e.g., Brekke et al., 2013) use GCMs that participated in CMIP3 and CMIP5 (note that there was not CMIP4). However, despite the inevitable improvements expected in GCM performance with time (Reichler and Kim, 2008) it turns out that the CMIP3 and CMIP5 projections are remarkably similar after accounting for the differences in GHG emissions scenarios (Knutti and Sedláček, 2013).
Using this climate model framework, we evaluate the hydrologic (watershed discharge, surface runoff, snow melt/accumulation), water quality (NO 3 -, dissolved P, sediment-bound P, total P, sediment export/yield), and N 2 , and N 2 O emissions implications of a changing climate. Results are presented at the aggregate watershed level for each individual climate model on an annual basis, for the ensemble model mean (each model run through SWAT-VSA and averaged) and range (minimum and maximum result derived from each climate model run through SWAT-VSA), and on a seasonal basis). Results are further broken down by developing quantile metrics to determine deviations from historic values. To do this we employ Q 10 and Q 90 metrics for all watershed outlet responses analyzed (e.g., flow, NO 3 -, dissolved P, total P, sediment export/yield). The Q 10 metric is the response value exceeded 90% of the time; a measure of low flows, and the Q 90 metric is the response value exceeded 10% of the time, a measure of high flows. We then assess the responses of agricultural HRUs, or fields to better explain the finer scale responses observed at the watershed level.

Selection of extreme wet/dry events
Extreme weather anomalies (wet and dry periods) from the biascorrected NARCCAP data were selected to test the effect of changes in soil moisture on nutrient cycling and the production of N 2 and N 2 O. Years with wet and dry periods were selected by developing indices of climate anomalies (Karl et al., 1999).To select years with dry periods, the maximum length of a dry spell, defined as the maximum number of consecutive days where daily precipitation is b1 mm, was counted for each year of both the historical and future climate periods. The average number of consecutive dry days across all climate models was 14 days with a standard deviation of 4 days.
Wet events were defined as the maximum number of consecutive days with daily precipitation depths N 1 mm, these were counted for both the historical and future climate periods. However, since there were a high number of wet events (according to our definition above) compared to the number of dry events, we reselected the top three wettest events in the time series. Model output corresponding to these events were extracted from the model runs analyzed for the change in the rate of N 2 and N 2 O production, and soil moisture relative to the historical climate data. This analysis focuses on field scale responses and on specific times (and not representative of the long term responses discussed above).

Watershed model performance assessment
According to performance criteria recommended by Daniel et al. (2015) and Daniel et al. (2007), the results indicate that the model predicted the observed flow, NO 3 -, and DP export well in both the calibration and evaluation periods (Table 1), with a low PBIAS and good NSE performance metrics. The model slightly under predicted streamflow, NO 3 -, and sediment and over predicted DP during the calibration period while slightly over predicted flow and NO 3 -and under predicted DP during the evaluation period (Table 1). Sediment predictions displayed slightly more bias during the evaluation period (note sediment was uncalibrated). N 2 O emissions were over predicted during the evaluation period, but the performance still falls within the acceptable range proposed by Daniel et al. (2015) particularly for daily model results.

Summary of future climatic projections
Precipitation projections from the seven climate models ranged from slight decreases in mean annual precipitation for two of the models (CRCM-CCSM and WRFG-CGCM3 models) to increases predicted by the remaining models (Table 2). Averaged across the seven models, mean annual precipitation increases by 3.7% compared to the historical period. All climate models predict an increase in precipitation intensity with the exception of HRM3-GFDL model, with the model mean intensity increasing by 5.1%. All models also predicted an increase in both the maximum and minimum temperature of 2.5°C. This increase in air temperature results in an average increase in the soil temperature of 1.9°C.

Climate change effects on hydrology
Climate change is predicted to cause small changes in the magnitude, and substantial changes in the timing of streamflow in the WE38 watershed. Ensemble model results suggest a small decrease of 0.7% in mean annual flow (Fig. 2a), but there is considerable variability among individual climate model predictions (−16.1% to +34.2%). The winter and spring flows (Fig. 3a) generally increases due to the increase in surface runoff (Fig. 4a) and the reduction in the ratio of snow to rain during the winter (Fig. 4b). Flow reductions in late spring and early summer (Fig. 3a) are primarily due to greater evapotranspiration (Fig. 4c).
The CRCM-CGCM3 and ECP2-GFDL models predict the greatest increase in annual flow (12.9% and 34.2%, respectively, Fig. 2a), while the CRCM-CCSM, HRM3-GFDL, RCM3-GFDL, MM5I-CCSM, and WRFG-CGCM3 models predict the greatest decreases in annual flow (1.4 to 16.1%, Fig. 2a). However, these changes appear to be driven by changes at extreme flows (e.g., all of seven models suggest decreases or no changes in Q 10 flows, and increases or no change in Q 90 flows, Fig. 2a).

Climate change effects on water quality
The ensemble model mean predicts an increase in mean annual NO 3 -watershed level export of 6.5% (Fig. 2b). Similar to streamflow, there is considerable variability among annual model predictions (−0.1% to +18.9%). All models suggest a substantial increase in Q 90 NO 3 -export and a decrease in Q 10 NO 3 -export (Fig. 2b), indicating substantially more variability than the historic record displayed. Peak NO 3export timing does not change considerably, historically occurring in March, during snowmelt, late May, following fertilizer application and again in December following crop harvest (Fig. 3b), but winter export increases significantly, from 60 to 80 kg d −1 to 85-105 kg d −1 . The historic summer/fall (June-November) NO 3 -export is substantially reduced under future conditions, primarily a result of decreased streamflow from reduced runoff (Fig. 4a) that result from drier summer soils, while the increase in winter/spring is due to greater runoff (Fig. 4a) and nitrification (Fig. 4g).
Model mean results indicate that mean annual DP export from the watershed will increase approximately by 16.6% (Fig. 2c). Similar to streamflow, there is considerable variability among the mean annual change predicted by the models for DP export (−1.7% to +55.9%). All models suggest substantial increases or no changes in Q 90 and decreases or no changes in Q 10 DP export (Fig. 2c). DP export increases across all winter/spring months, generally by 20-30% (Fig. 3c). Peak DP export continues to occur in April, but increase from 0.5 kg d −1 to 0.6 kg d −1 . Annual DP yield from agricultural fields increases 10% (Fig. 4f), despite a decrease in soil P mineralization (Fig. 4h), the conversion of insoluble organic P to soluble inorganic P. This is due to a significant increase in surface runoff (Fig. 4a) and sediment-bound P in the watershed (Fig. 4e).
The mean annual watershed level TP export (sum of DP, sedimentbound P, and organic P) is predicted to increase (10.8%), however, there is a wide range of predicted responses (−0.2% to +43.2%, Fig. 2d). All model predictions result in a decrease or no changes in Q 10 and increases in Q 90 TP export (Fig. 2d). Peak TP export, historically occurring in February-March, coinciding with the spring snowmelt flush, shifts to April-May, when fields in the watershed have been tilled (Fig. 3d). Peak TP export also increases from a mean of 2.0 kg d −1to 2.5 kg d −1 , although the maximum predicted by individual models was far greater, up to 7.0 kg d −1 (Fig. 3d).
At the agricultural field scale, almost all models agree that sedimentbound P yield will increase substantially, driven by large winter/spring increases (Fig. 4e). Mean annual increases range from 2.1% to 47.4% (Fig. 4e), with a model mean increase of 20.9%. Increases in the winter/spring range from 25.7% to 78.2% (Fig. 4e), while decreases are predicted during the summer/fall (−6.0%). The increase in peak export, particularly in the winter/spring is due to the combined effect of increased surface runoff (Fig. 4a) and rainfall intensity (Table 2), coupled with reduced snow cover (Fig. 4b).
Mean annual watershed level sediment export increases 6.7% during the future period, again, with considerable variability among models (−4.5 to +35.8%, Fig. 2e). There are also substantial decreases or no changes in the Q 10 (−71%) and increases in Q 90 (17.1%) sediment export (Fig. 2e). Much of the annual increase is due to increased sediment export in the winter/spring; although there is decrease in summer sediment export (Fig. 3e). Peak sediment export increases from 1.0 ton d −1 , historically occurring in late February-early March, to 1.5 tons d −1 in April-May, again, coinciding with tillage of the agricultural fields and reduced ground cover (Fig. 3e).
Mean annual sediment yield from agricultural fields in the watershed increases by 8.0% (Fig. 4d). All models predicted substantial increases in sediment yield during the winter/spring (6.6% to 43.5%) and most models predicted decreases in the summer/fall (four of seven models, Fig. 4d). This is consistent with and of similar magnitude as the changes in surface runoff predicted for both periods (Fig. 4a). Note that the increase in sediment-bound P (Fig. 4e) is proportional to the increase in sediment yield (Fig. 4d) during both periods.  a The change in the precipitation intensity is calculated by summing the precipitation in time series divided by the total number of days with precipitation in time series. b The change in soil temperature is taken directly from the SWAT model output.  Fig. 4i shows mean annual decreases in the emission of N 2 O (−9.9%), which is driven by significant decreases in summer N 2 O emissions (−27.0%), primarily due to lower soil moisture inhibiting microbial processes. Conversely, N 2 O emissions increase from December through late March due to increases nitrification (Fig. 4g) and greater soil moisture.

Emissions of N 2 O and N 2
N 2 emissions similarly decrease on an annual basis (9.4%, Fig. 4j), again, driven by significant decreases is during the summer/fall (−36.6%) due to drier soils. The slight increase in winter/spring N 2 emissions are driven by increased soil nitrification (Fig. 4g), which provides more labile NO 3 -(the reactant for denitrification) coupled with increased soil temperatures and moisture (Table 2), which speeds up microbial denitrification.  Historical (1975Historical ( -1998 and future (2045-2068) average-by-day from the ensemble model mean and the ensemble model range for streamflow (a), and the export of nitrate (b), dissolved P (c), total P (d), and sediment (e) for the seven regional climate models.
The greatest increase in N 2 O emissions occur in October-March (Fig. 5a), which coincides with a period of high watershed discharge (Fig. 3a), indicative of saturated soils. Similar to N 2 O, the ensemble model mean predicts large decreases in emissions of N 2 from late March to late November but then predicts increases from December through March (Fig. 5b). Decreased emissions of both N 2 O and N 2 during summer/fall are driven by higher summer soil temperatures (Table 2), and decreases in the nitrification rate in the soil (Fig. 4g), making less N available for denitrification (Figs. 4i and j). Increased soil temperatures (Table 2) similarly increase the microbial respiration rates, which fosters denitrification. Thus, despite the increase in microbial respiration rate, less N 2 O and N 2 emissions are predicted due to less available labile NO 3 -during summer/winter (Fig. 3b), which is the main driver of denitrification. Fig. 6 shows the spatial model predictions of N 2 O and N 2 emission under wet (Fig. 6a) and dry (Fig. 6b) conditions. The high emissions areas of both N 2 O and N 2 tend to coincide with areas of high soil moisture irrespective of whether the average watershed condition is wet or drier. That is, hotspots of N 2 and N 2 O emissions are generally found in areas of higher soil moisture. The management influence is also clear in Fig. 6, where high soil moisture areas do not always result in areas of high emissions; high emission areas are also found on landuses that receive substantial fertilizer or manure additions, such as corn, increasing the N available to drive denitrification.

Discussion
Ensemble modeling highlights how climate change and climate anomalies are likely to affect water availability and water quality in the Chesapeake Bay region and the level of uncertainty in those predictions. The results from this analysis indicate the region is likely to experience a slight decrease in mean annual flow although the magnitude of the change differs among the climate models (Fig. 2a). However, an increase in flow is predicted during the winter/spring (10.6%) while a significant decrease in flow is predicted during the summer/fall (−29.1%, Figs. 2a and 3a). These results corroborate findings by others in the region. For instance, Hayhoe et al. (2007), Rob et al. (2000), Lu et al. (2015), and Najjar et al. (2010) all predicted increases in seasonal Fig. 4. Annual and seasonal changes (%) for surface runoff (a), snowmelt (b), actual evapotranspiration (c), sediment yield (d), sediment-bound P export (e), dissolved P export (f), nitrification rate (g), phosphorus mineralization rate (h), N 2 O emission rate (i), and N 2 (j) emission rate averaged across agricultural HRUs for seven climate models relative to historical scenario  for the future time period (2045-2068). streamflow variability with decreased flow during the summer due to increased evapotranspiration (driven by higher temperatures) and increased flow in the winter and spring due to increased winter precipitation and reduced snowpack.
The predicted decline in flow during summer from watersheds such as WE-38 could affect the ecology of the Chesapeake Bay estuary; in particular reducing dissolved oxygen, instream habitat, and increasing the potential for harmful algae blooms (Meyer et al., 1999;Neff et al., 2000). The increase in flow during winter/spring occurring on fields that have been harvested and thus have less plant cover, could result in increased sediment delivery to both reservoirs in the watershed and to the estuary (Cerco, 2016). Furthermore, the increase in winter flow is indicative of higher soil temperatures (Table 2), and reduced snowpack (Fig. 4b) in the watershed. These changes could result in producers in the regions shifting planting, field tillage, and fertilizer operations to earlier in the year to take advantage of warmer temperatures (Sims and Coale, 2002), which could change the timing of nutrient and sediment export from the watersheds in the region.
Sediment export tends to be a function of streamflow; mean annual sediment export increased across most models (Fig. 4d). This increase in sediment export from landscape is due to increased runoff and precipitation intensity during winter and spring when there is reduced groundcover following harvest (Fig. 4d). The increased future sediment export is primarily derived from landscape sources (as opposed to channel sources). This is determined by subtraction: the HRU sediment yield was subtracted from the sediment yield at the outlet of the watershed. The ratio of channel-originated sediment to the HRU-originated sediment is 0.18 for the historic period and 0.20 for the future period. However, over time increased deposition of landscape sediment in channels could increase the contribution of the channel to overall export. Increased sediment deposition in the channels could also alter stream geomorphology, increase flooding and stream water temperature (Rehg et al., 2005), reduce light penetration, and affect submerged aquatic vegetation (Bilotta and Brazier, 2008).
Phosphorus export exhibited substantial variability across models; with some climate models predicting slight increases in mean annual P export while others predict decreases (Figs. 2c and d and 4e and f). The ensemble model mean predicted an increase in mean annual P export; which is counter to predictions made by Roberts et al. (2009) for the Chesapeake Bay, who predicted a 19% reduction in P using the Spatially Referenced Regression On Watershed Attributes (SPARROW) model, but similar to results by Chang (2004) for southeastern Pennsylvania, who predicted increased P export. Our results suggest that the increase is largely driven by increases in sediment-bound P export during the winter/spring (Fig. 4e), and increased precipitation intensity (Table 2) and not from increased soil P mineralization (Fig. 4h, P mineralization actually decreases). The increase in P export could result in eutrophication and depletion of dissolved oxygen in the stream (Boesch et al., 2001). Furthermore, the increase in spring DP and TP export is particularly troubling for the estuary as P delivered in the spring flush tends to set up more aggressive summer eutrophication conditions (Boesch et al., 2001).
Most models used in this analysis predicted an increase in mean annual and winter/spring NO 3 -export, due to significant increases during the winter/spring (+17.3%) but significant decreases during the summer/fall season (−29.1%), similar to results by Chang (2004) for southeastern Pennsylvania. The annual increase in NO 3 -export is due to substantial increases in nitrification during the winter/spring (Fig. 4g), and increased runoff/soil moisture. Increased nitrification creates a greater pool of N available for transport as NO 3 -is more mobile in the soil than organic N or NH 4 -(the reactants in N mineralization and nitrification, respectively). During the winter/spring, the increase in nitrification ( Fig. 4g) provides more NO 3 -to the system than can be utilized by denitrification because this system is C limited (b5:1C:N ratio). This process is driven by increased temperatures, soil moisture levels, and soil N levels (Agehara and Warncke, 2005;Guntiñas et al., 2012). In terms of spatial variability, these results corroborate the hot moment/ hot spot theory of N emissions. Indeed, Fig. 6 shows those areas of the landscape with greater soil moisture levels are the primary sources of elevated N 2 O emissions. However, the spatial heterogeneity also depends on soil NO 3 -and C availability, soil temperature, and management practices (Wagena et al., 2017).
These results indicate that agricultural soils continue to serve as a source of NO 3 -export, even without fertilization, for a long period of time (average soil total N is 16,667 kg/ha) posing a threat to water quality. However, there are agricultural conservation practices that could reduce NO 3 -export from the watershed. For example, Noellsch et al. (2009) proposed slow-or controlled-release N fertilizers to reduce NO 3 -leaching losses, Carpenter et al. (1998) also recommended landscape management (e.g., riparian vegetation, restoration of wetlands, and floodplains) to increase denitrification that decreases NO 3 -export. Alternatively, agriculture in the region may adapt to the new conditions, for instance by planting earlier to take advantage of increased soil temperatures and soil moisture in the spring, or by introducing new crop varieties that may be adapted to the higher temperatures, altered soil moisture levels, and/or nutrient cycling rates.

Conclusions
Climate change will likely alter water availability and quality. Results indicate that the future climate projected by most of the climate models used in this study will result in increased flow, sediment, P, and NO 3 -export during the winter and spring, while decreases are expected during the summer and fall, although there is considerable variability among models. The increase in NO 3 -export due to nitrification and N mineralization could potentially threaten water quality. Increased P export is driven primarily by increased sediment-bound P export from the watershed, driven by increased surface runoff and rainfall intensity coinciding with periods of reduced soil plant cover. These results suggest biogeochemical cycles could be substantially altered, which may ultimately result in changes to agricultural management practices. The results of this effort will also help managers to develop strategies to maintain crop productivity, and water quality. Sustainable intensification should focus on the timing of planting and fertilizer applications so that crops are actively growing when fertilizers are applied while reducing nutrient export. Precision agriculture and critical source area models can be utilized to determine where and how much fertilizer should be applied. With proper planning, agricultural management could increase crop yields, reduce nutrient exports, and improve water quality.