Continuous separation of land use and climate effects on the past and future water 1 balance

19 Understanding the combined and separate effects of climate and land use change on the water cycle is 20 necessary to mitigate negative impacts. However, existing methodologies typically divide data into 21 discrete (before and after) periods, implicitly representing climate and land use as step changes when in 22 reality these changes are often gradual. Here, we introduce a new regression-based methodological 23 framework designed to separate climate and land use effects on any hydrological flux of interest 24 continuously through time, and estimate uncertainty in the contribution of these two drivers. We present 25 two applications in the Yahara River watershed (Wisconsin, USA) demonstrating how our approach can 26 be used to understand synergistic or antagonistic relationships between land use and climate in either the 27 past or the future: (1) historical streamflow, baseflow, and quickflow in an urbanizing subwatershed; and 28 (2) simulated future evapotranspiration, drainage, and direct runoff from a suite of contrasting climate and 29 land use scenarios for the entire watershed. In the historical analysis, we show that ~60% of recent 30 streamflow changes can be attributed to climate, with approximately equal contributions from quickflow 31 and baseflow. However, our continuous method reveals that baseflow is significantly increasing through 32 time, primarily due to land use change and potentially influenced by long-term increases in groundwater 33 storage. In the simulation of future changes, we show that all components of the future water balance will 34 respond more strongly to changes in climate than land use, with the largest potential land use effects on 35 drainage. These results indicate that diverse land use change trajectories may counteract each other while 36 the effects of climate are more homogeneous at watershed scales. Therefore, management opportunities to 37 counteract climate change effects will likely be more effective at smaller spatial scales, where land use 38 trajectories are unidirectional. 39


Introduction
Climate and land use (which we define broadly to include land use, land cover, and land management) are two major drivers of global hydrological change (Foley, 2005;Steffen et al., 2015;Vörösmarty et al., 2000).While economic, governmental, and social pressures may be exogenous to a watershed, land use can be controlled by decision-making at local levels (individual, city, county, and state).In contrast, climate change is driven by global emissions and requires a coordinated effort well beyond an individual watershed to address.Therefore, land use decisions may be a viable path to mitigating undesirable impacts of climate change on the water cycle at watershed scales.While several point-based studies have found significant impacts of land use change on the water balance (Giménez et al., 2016;Nosetto et al., 2012;Scanlon et al., 2005;Twine et al., 2004), most watershedscale studies have found that the impact of climate change on hydrology outweighs that of land use change, particularly where there are strong changes in precipitation (Chawla & Mujumdar, 2015;Jiang et al., 2015;Li et al., 2009;Mango et al., 2011;Pumo et al., 2017;Tao et al., 2014;Wu et al., 2015;Yang et al., 2017).Thus, there is a growing acknowledgment that the impacts of land use change are superimposed on a larger climate trend, and can either amplify or partially counteract the impacts of climate change (Gyawali et al., 2015;Juckem et al., 2008;Martin et al., 2017;Shi et al., 2012;Tomer & Schilling, 2009;Wang & Stephenson, 2018;Zhang et al., 2016).In particular, land use may be most important locally (Frans et al., 2013;Haddeland et al., 2007;Peterson et al., 2011;Xu et al., 2013), during wet/dry extremes (Villarini & Strong, 2014), or where significant infrastructure projects (e.g.dams) occur (Wu et al., 2012;Ye et al., 2003).
However, existing statistical approaches to separate the effects of climate and land use change have three major limitations.First, existing statistical methodologies often implicitly treat land use and climate effects as step-changes by dividing datasets into two or more discrete time periods (e.g."before" and "after") (Gupta et al., 2015;Li et al., 2009;Tomer & Schilling, 2009;Wang & Hejazi, 2011;Wang & Stephenson, 2018;Xu et al., 2013;Zhang et al., 2016).This assumption may be problematic because land use and climate often change continuously and in tandem, with gradual hydrological impacts (Jiang et al., 2015;Marhaento et al., 2017).Therefore, there is a need for improved methods to separate the impacts of climate and land use and their interactions in a continuous time series of hydrological data.
Second, with some exceptions (e.g,.Tao et al., 2014), previous studies have primarily focused on disentangling the relative importance of climate and land use on historical streamflow data, particularly Budyko-based methods which assume that ET is equal to the difference between precipitation and runoff (Jiang et al., 2015;Renner et al., 2012;Tang & Wang, 2017;Wang & Hejazi, 2011;Wang & Stephenson, 2018;Xu et al., 2013).In order to adequately understand and address the impacts of climate and land use on water resources, tools are needed to quantify the impacts of these drivers on the complete water cycle which includes hydrologic fluxes of evapotranspiration [ET], drainage, and runoff.
Third, most existing statistical regression-based approaches rely on multiple linear regression (MLR) with precipitation and potential ET as predictors (e.g.Ahn & Merwade, 2014;Huo et al., 2008;Jiang et al., 2011;Ye et al., 2003;Zhang et al., 2016).This a priori assumption about relevant variables to hydrological fluxes ignores other potentially important meteorological drivers, for example vapor pressure deficit, incoming solar radiation, and wind speed (Vicente-Serrano et al., 2014;Zhou et al., 2014;Zipper et al., 2017b).Additionally, MLR assumes no multicollinearity between input variables, but in reality meteorological variables may be correlatedfor instance, wetter months are often cooler.Thus, there is a need for approaches which both consider a greater diversity of potential input variables for regression, while simultaneously addressing multicollinearity issues that compromise statistical predictability.Approaches to transform input data to maximize variance and orthogonality include (1) principal components regression (PCR), which maximizes variance of predictor variables; and (2) partial least squares regression (PLSR), which maximizes the variance of the predictor variables relative to the response variable (Hadi & Ling, 1998;Hwang & Nettleton, 2003;Wang, 2012;Wolter et al., 2008).However, these regression models have not been rigorously tested or intercompared for hydrologic applications.
To meet this challenge, we introduce a new methodological framework to quantify the impacts of climate and land use change on any measured or modeled hydrological flux continuously through time.Using this approach, we demonstrate that PLSR outperforms both PCR and the more commonly used MLR (differences between these regression models are described in Section 2.1.3).We then apply this method to answer the question, to what degree can land use amplify or counteract climate-induced changes to the water balance of a watershed?Using both historical data and simulated results from diverse future scenarios for streamflow, ET, drainage, and direct runoff, we provide insight into the degree to which land use can be used as a local tool to maintain a watershed within a desired hydrological operating space in the context of an uncertain future climate.

Statistical model description
In brief, our new method develops statistical relationships between meteorological variables and any hydrological flux (HF) of interest during a baseline period.These statistical relationships are then applied to climate data outside of the baseline period, which we refer to as the prediction period.Predictions are used to estimate the HF changes resulting from climate change relative to the baseline period at a monthly resolution, and residuals from predictions are attributed to land use.A general overview of the method is presented in Figure 1.The specific type of regression model used within our methodological framework may vary; in this study, we compare partial least squares regression (PLSR), principal components regression (PCR), and multiple linear regression (MLR) models.While the analysis in this study is done at a monthly timestep consistent with other regression-based studies separating climate and land use effects (Ahn & Merwade, 2014;Schottler et al., 2014;Xu et al., 2013;Ye et al., 2003;Zhang et al., 2016), the method may be applied at other time resolutions as long as reliable input data and regression relationships can be developed.
We describe this method for a generic HF in the present section (2.1), and then separately apply it to historical streamflow data (Section 2.2.2); and simulated future ET, drainage, and direct runoff (Section 2.2.3).

Generating baseline relationships
To estimate the relative contributions of climate and land use to change in a given HF, a baseline period must be identified from which relative changes are then calculated.This baseline period should represent a period of time in which land use is relatively static, so that variability in the HF is driven primarily by meteorological processes.First, we use significance pruning to select predictor variables for the HF of interest from a suite of candidate predictor variables within the baseline period.It is important that candidate predictor variables are (a) available over the entire period of interest; and (b) controlled by climate, not land use.In our applications (Section 2.2), candidate variables include a variety of measured and derived meteorological variables (e.g.precipitation, temperature, reference ET).For each month, candidate predictor variables were mean-centered and scaled to a unit standard deviation.We retain the subset of candidate predictor variables that have a significant linear relationship with the HF (a significance threshold of p<0.10 was used to err on the side of variable retention).Importantly, this approach means that the retained predictor variables are allowed to differ by each month and HF; for example, incoming solar radiation may be a more important predictor for ET than direct runoff.
To avoid potential overfitting and eliminate collinearity between predictor variables, we use PLSR to transform the retained variables to factors which we use to predict the HF of interest.To determine the factors used for regression models, we use the factors explaining a cumulative 80% of total variance in the scaled predictor variables.The selected factors are used as input to a multiple linear regression equation of the form: HFPLSR,m,y = C0 + C1*F1,m,y + C2*F2,m,y + … + Cn,m*Fn,m,y + ε {Eq.1} where HFPLSR,m,y is the hydrological flux for month m and year y predicted by PLSR, Cx are regression coefficients, Fx,m,y are factors, and ε is an error term assumed to be normally distributed and centered on 0. We use a permutation-based, split-sample approach to estimate model fit and uncertainty, in which we run the PLSR 250 times randomly sampling 75% of the baseline period for model calibration while retaining 25% for model validation (Zipper & Loheide, 2014).This approach provides 250 unique sets of regression coefficients for each month and HF.

Calculating climate and land use contributions to change
The statistical relationships for the baseline period are then applied to the rest of the hydrological time series (the prediction period).While both of our applications of the method (Section 2.2) have a baseline period before the prediction period, this method can also use modern conditions as the baseline period and apply statistical relationships into the past to quantify the relative contribution of historical land use and climate change to a given HF, as long as land use in the baseline period is relatively stable.Using the permutation-based approach described above, we have 250 estimated values of each HF for each year and month within the prediction period.
To separate the relative contribution of climate and land use, we adopt the common assumption that these two factors can explain all variability in a given HF relative to the baseline period: climate encompasses all changes due to drivers exogenous to the study system (in our case, the watershed), and land use encompasses all changes due to drivers endogenous to the study system (Ahn & Merwade, 2014;Duan et al., 2017;Gao et al., 2016;Huo et al., 2008;Jiang et al., 2011).Therefore, changes in land management such as irrigation or fertilization practices are included in the land use category, while longer-term climatic oscillations (e.g.PDO, ENSO) would be included in the climate category.
For each HF, we calculate the total change relative to the baseline period (ΔHFTotal,m,y) as: ΔHFTotal,m,y = HFm,y -HFm,baseline {Eq.2} where HFm,y is the measured or modeled hydrological flux for month m and year y, and HFm,baseline is the mean HF for that month during the baseline period.The total climate contribution to change (ΔHFClimate,m,y) can then be expressed as: ΔHFClimate,m,y = HFPLSR,m,y -HFm,baseline {Eq.3} where HFPLSR,m,y is the PLSR-estimated value for month m and year y.Finally, the land use component of change (ΔHFLU,m,y) is calculated as: ΔHFLU,m,y = ΔHFTotal,m,y -ΔHFClimate,m,y = HFm,y -HFPLSR,m,y {Eq.4} Note that any of the ΔHF variables can be positive or negative, corresponding to an increase/decrease in that HF relative to the baseline period.This framework allows us to quantify not just the overall change relative to the baseline period for each month, but also under what conditions the effects of land use and climate are antagonistic (ΔHFLU,m,y and ΔHFClimate,m,y have opposite signs) and under what conditions the effects of land use and climate are synergistic (ΔHFLU,m,y and ΔHFClimate,m,y have the same sign).

Comparing among regression models
While the overall methodological framework we introduce in Section 2.1.1 and 2.1.2uses PLSR to generate statistical relationships, our approach will work with any regression model capable of accurately predicting HF at a monthly resolution.To determine the importance of the choice of regression model, we repeated our analysis using principal components regression (PCR) and multiple linear regression (MLR) approaches.We conducted this comparison using historical discharge data from the Pheasant Branch subwatershed, which is described in Section 2.2.2.
When our methodological framework was adapted to use PCR and MLR regression models, we used the same general approach described in Section 2.1.1 and 2.1.2:statistical relationships were generated between meteorological variables and a response variable for a baseline period, then applied monthly through time for a prediction period.Total change was calculated as the difference between a monthly flux and the mean of the baseline period (e.g.Eq. 2), changes due to climate were calculated as the differences between the mean of the baseline period and each statistically-predicted monthly value (e.g.Eq. 3, replacing HFPLSR with HFMLR or HFPCR for the MLR and PCR approaches, respectively), and changes due to land use were calculated as the difference between total change and changes due to climate change (e.g.Eq. 4, replacing HFPLSR with HFMLR or HFPCR).
All three regression models used the same significance pruning approach to select candidate predictor variables (Section 2.1.1),and therefore these were the same for all methods.To determine the principal components (PCs) used for PCR, we selected the PCs explaining a cumulative 80% of total variance in the scaled predictor variables, as well as any other PC which has both a significant linear relationship with the HF (p<0.10 as above), and explains >1% of total variance in input variables (to avoid spurious correlations), as low-ranked PCs can be important for PCR (Hadi & Ling, 1998;Jolliffe, 1982).For MLR, we used the n best predictor variables (as measured by linear R 2 ), where n is equal to the number of PCs retained in PCR.For all methods, we calculated model fit using 250 permutations where data were split into 75% calibration/25% validation samples as described in Section 2.1.1.

Study area
We applied the approach described in Section 2.1 to the Yahara River Watershed (YW; area=1344 km 2 ), Wisconsin, USA (Figure 2).The YW is an urbanizing agricultural watershed, and thus is a useful analogue for human-influenced watersheds throughout the US Midwest and the world (Carpenter et al., 2015b).The water resources of the YW are stressed by various land use and climatic drivers of change including (1) an expanding urban core (Madison, the state capital), leading to changes to the water and energy balance (Schatz & Kucharik, 2014;Zipper et al., 2016Zipper et al., , 2017b)); (2) widespread fertilized row-crop and dairy agriculture contributing to erosion and nutrient loading (Carpenter et al., 2015a;Lathrop & Carpenter, 2013;Motew et al., 2017;Qiu & Turner, 2013, 2015); and (3) a long-term trend of increasing precipitation with more frequent extreme precipitation events in recent decades, leading to both groundwater and surface water issues (Figure S1; Booth et al., 2016a;Gillon et al., 2016;Usinowicz et al., 2017).Due to these stresses on the water cycle, improving the understanding and management of climate and land use effects on water resources is a key goal cutting across hydrological, ecological, and social research in the YW (Gillon et al., 2016;Motew et al., 2017;Qiu et al., 2017Qiu et al., , 2018;;Wardropper et al., 2015).We separately investigated the YW's past (Section 2.2.2) and future (Section 2.2.3) to quantify how the water cycle of the YW has changed historically and may continue to change under a variety of future scenarios.
For clarity, throughout the paper the term "discharge" is used to refer to total streamflow as measured at a gauging station converted to units of depth after dividing by total watershed area; discharge can be separated into "quickflow" and "baseflow" components (Schwartz & Smith, 2014)."Direct runoff" is used to refer to overland flow calculated at the grid cell resolution from Agro-IBIS output (Section 2.2.3).

Historical discharge analysis
For historical analysis, we focused on the Pheasant Branch Subwatershed (PBS; 44.24 km 2 ; Figure 2) which drains the northwest portion of the YW including parts of the municipalities of Madison and Middleton.We selected the PBS for detailed analysis because it has the longest available record of discharge data (1974-present) in the YW upstream of the chain of lakes (Figure 2), which buffer the impacts of climate on streamflow at the monthly scale of analysis used here.Since 1974, there have been well-documented changes in land use (urbanization, including the connection of former internally-drained basins to the streamflow network), water governance (stringent infiltration requirements for new developments), climate (increased precipitation), and flood peaks (increasing discharge) (Gebert et al., 2012).
We applied the our method using monthly discharge data from the USGS National Water Information Service gauging station 05427948 (U.S. Geological Survey, 2017) for the period July 1974-December 2016 (a total of 42 years and 6 months).We defined the baseline period as the first half of the available streamflow data (July 1974-December 1995; 21 years and 6 months), and the prediction period as the second half of available discharge data (January 1996-December 2016;21 years).This breakpoint also roughly corresponds with an observed shift in historical streamflow beginning in 1993, which has been attributed to increasing precipitation and urbanization within the PBS (Gebert et al., 2012).
Predictor variables for the PBS were either measured or derived from the Madison Airport Global Historical Climatology Network-Daily (GHCN-D) site (USW00014837; 43.14°N, -89.35°E) (Menne et al., 2012).Directly measured variables were daily precipitation, maximum temperature, and minimum temperature.Wind speed data was available for only part of the period of interest, and therefore we used mean values for a given day of year for the entire period.We estimated vapor pressure as the saturation vapor pressure at minimum daily temperature following Allen et al. (1998).We estimated daily incoming solar radiation using the Bristow-Campbell equation (Bristow & Campbell, 1984), which scales the topof-atmosphere solar radiation using an estimated transmissivity based on daily maximum and minimum temperature, as implemented in the EcoHydRology R package (Fuka et al., 2014).The Bristow-Campbell equation was calibrated to site conditions using observed incoming shortwave radiation data from the nearby Arlington Agricultural Research Station (43.31°N, -89.38°E; http://agwx.soils.wisc.edu/uwex_agwx/awon)for the period 1986-2016 (Figure S2).We calculated daily Penman-Monteith reference ET following the UN Food and Agriculture Organization method (Allen et al., 1998), and precipitation deficit as precipitationreference ET.
We then aggregated daily variables to a monthly set of candidate predictor variables: cumulative monthly precipitation, reference ET, and precipitation deficit [mm mo -1 ]; and mean daily minimum and maximum temperature [°C], incoming shortwave solar radiation [W m -2 ], wind speed [m s -1 ], relative humidity [%], actual vapor pressure, saturation vapor pressure, and vapor pressure deficit [kPa].Candidate predictor variables included both the month of interest, as well as the month of interest plus the preceding 1, 2, 3, 6, and 12 months by summing (cumulative variables) or averaging (mean daily variables).We also included monthly metrics associated with precipitation intensity, including maximum daily precipitation [mm], total days with precipitation, and total days with precipitation exceeding 12. 7, 25.4, 50.8, and 76.2 mm (0.5", 1", 2", 3") which are becoming more frequently exceeded in the YW recent decades (Gillon et al., 2016); and metrics allowing for nonlinear responses to precipitation, including squared monthly precipitation [mm], squared monthly precipitation deficit, and squared cumulative precipitation deficit for all lags.In total, there were 79 candidate predictor variables evaluated for each month.The retained variables for each flux are shown in Figure S3.
We also performed a parallel set of analyses for the quickflow and baseflow components of discharge in the PBS separated using a recursive digital filter (Eckhardt, 2005) within the Web-based Hydrography Analysis Tool (WHAT; Lim et al., 2005).All other analyses were repeated as described above using GHCN-D data from the Madison Airport.These results are presented in the Supplementary Information.
To supplement our permutation-based approach to estimating model uncertainty, we conducted two additional analysis for historical discharge in the PBS to quantify the sensitivity of our method to (i) the selection of a baseline period; and (ii) uncertainty in meteorological input data.For (i), we repeated all analyses while varying the end of the baseline period from 1992 to 1998 (1995 +/-3 years).For (ii), we repeated our analyses with a baseline period ending in 1995 using meteorological data from the Arlington Research Farm GHCN-D site (USC00470308), which is in the northern part of the Yahara River Watershed.Arlington was missing precipitation data for 478 days during the study period (a total of 15,525 days), which were gap-filled using precipitation from the Madison Airport GHCN-D site.

Biophysical model description
To investigate the extent to which climate and land use may impact different components of the water balance, we simulated a variety of plausible future scenarios for the YW using Agro-IBIS, a gridded, physically-based dynamic vegetation model including agroecosystems.Agro-IBIS simulates the complete carbon, energy, and water cycles (Foley et al., 1996;Kucharik et al., 2000;Kucharik, 2003;Kucharik & Brye, 2003).Recent updates to Agro-IBIS replaced the soil physics with those of HYDRUS-1D (Šimůnek et al., 2013), so that the soil water balance is solved using the pressure head-based form of the Richards' Equation (Soylu et al., 2014); and added erosion and phosphorus cycling, along with a suite of new land cover types, for the simulation of the YW (Motew et al., 2017).
In this study, we used the version of Agro-IBIS described in Motew et al. (2017), which simulates the YW at 220-m×220-m spatial resolution.Motew et al. (2017) calibrated and validated the hydrologic performance of the model via comparison with long-term streamflow records from six USGS gauging stations within the YW.Sediment/phosphorus transport and soil phosphorus concentrations were also validated against measurements.Previous validations of Agro-IBIS in the YW include comparisons against plot-scale measurements of soil moisture, soil temperature, leaf area index, aboveground net primary productivity, drainage, nitrogen cycling, and corn yield (Kucharik & Brye, 2003;Soylu et al., 2014;Zipper et al., 2015).In the interest of space, the reader is referred to the publications referenced above for additional information on the structure and validation of Agro-IBIS for the YW, though a validation with observations for the PBS is shown in the SI.

Climate and land use scenarios
Four scenarios, each with a unique climate and land use pathway, were developed to explore alternative social-political options for human action and socio-economic development in the YW for 2014-2070.Details of the storylines and biophysical drivers are presented in Carpenter et al. (2015b), Wardropper et al. (2016), and Booth et al. (2016b).The use of stakeholder-driven qualitative scenarios acknowledges the many potential paths climate and land use may take in the future, rather than focusing on a single forecasted future, and allows us to explore the degree to which climate and land use may interact under a variety of futures (Blöschl & Montanari, 2010).
Each of these four scenarios contains a separate set of land use and climate input data (Figure 3), as well as differences in the crop response to water stress representing agricultural biotechnology improvements.Full narratives, videos, and other information regarding the scenarios are provided in the abovereferenced publications and at yahara2070.org.A brief summary of key land use and climate drivers for each of the four scenarios follows: Accelerated Innovation (AI): AI explores a future in which technology is prioritized as a solution to climate change.Land use is characterized by expanding urban areas, with a relatively constant agricultural footprint.Climate change is the least extreme in this scenario, with warming of ~2°C by 2070 and more frequent heavy rainfall events.
Abandonment and Renewal (AR): AR explores a future in which society is unprepared for climate change.A mass exodus from the YW leads to a reduction in urban and agricultural land use, and the landscape primarily returns to natural vegetation.Climate change is the most extreme in this scenario, with warming of 5.5°C by 2070 and a period of extreme heat waves and floods in the 2030s.
Connected Communities (CC): CC explores a future in which sustainability and community become global priorities.Urban land use stays relatively constant, but agricultural land shifts away from row-crop agriculture to pasture and crops used directly as food (e.g.vegetables and small grains).Climate change in this scenario is intermediate between AI and AR, with 3.5°C warming by 2070 and both heavy rainfall and drought becoming more common.
Nested Watersheds (NW): NW explores a future in which governance is focused around national-scale water security.Urban land use remains relatively constant, but row-crop agriculture decreases as natural ecosystems are prioritized for water quality protection.Climate in this scenario is comparable to CC, with 4°C warming by 2070 and more frequent precipitation extremes.We used the approach described in Section 2.1 to evaluate climate and land use impacts on three HFs: ET, drainage, and direct runoff.These variables were averaged monthly over all terrestrial grid cells in the YW based on simulation output from the calibrated Agro-IBIS model of the YW (Motew et al., 2017).
The model was spun-up for 200 years using randomly selected meteorological years from the 1986-2013 period to equilibrate water, energy, carbon, nitrogen, and phosphorus cycles.The 1986-2013 period, during which land use and climate were the same for all scenarios, was used as the baseline period.
For the prediction period (2014-2070), we simulated a factorial combination of all land use and climate scenarios (16 total simulations) in order to provide a wide range of scenarios to evaluate interactions between land use and climate change.We used the same meteorological predictor variables as in the historical streamflow analysis (Section 2.2.2), though in this case they were watershed averages derived from spatially variable gridded meteorological input datasets (Booth et al., 2016b).The retained variables for each HF are shown in Figure S4.As in the historical streamflow analysis, we fit the PLSR model using 250 randomly sampled permutations of calibration/validation data which divide the baseline period into 75%/25% of available years.
For direct comparison with analysis of the PBS, we also extracted modeled monthly direct runoff for the 1974-2016 period from all grid cells within the PBS.For the Agro-IBIS spin-up, which includes 1974-1985, spatially distributed precipitation data were not available so randomly sampled meteorological years from the period 1986-2013 were used.For the 2014-2016 period, climate and land use from the AI scenario were used, though all scenarios are similar during this period.Therefore, we used the 1986-2013 period to compare Agro-IBIS' direct runoff performance for the PBS with quickflow derived from baseflow separation, and the entire 1974-2016 period for separation of climate and land use effects (with a 1974-1995 baseline period, as in the historical analysis).

Comparison among regression models and validation
While all of the regression models perform acceptably (Nash-Sutcliffe Efficiency [NSE] > 0; Moriasi et al., 2007) at both monthly and annual timescales, validation performance of PLSR is superior to both PCR and MLR (Figure 4).Notably, PLSR better predicts streamflow peaks (e.g.1993), leading to a monthly NSE of 0.716 (compared to monthly NSE < 0.5 for both PCR and MLR) and annual NSE of 0.878 (compared to annual NSE < 0.7 for PCR and MLR).Furthermore, the PLSR method is the only statistical method to have a root mean squared error (RMSE) <10% of the observed range of variability at either the monthly or annual resolution.PLSR also has a smaller range of predictions during both the validation and prediction period than PCR and MLR.However, even PLSR has a tendency to underestimate high monthly values of discharge (e.g. in the late 1970s), though to a lesser degree than either PCR or MLR.Based on the results of this comparison, we elected to use PLSR for the remainder of the analysis presented in the manuscript.

Climate and land use impacts on discharge
Within the baseline period , the method forces changes in discharge due to overall, climate, and land use effects to a mean of 0 mm, as the baseline period is the datum from which changes are calculated during the prediction period.During the baseline period, there is a slight but not significant trend in overall changes in discharge and climate-induced changes in discharge of 1.6 mm yr -1 (p=0.2), and no trend in land use-induced changes (slope=0 mm yr -1 ).This lack of a land use trend during the baseline period indicates that there is no trend in the residual of the PLSR relationships, lending support to our baseline period selection.
Within the prediction period , there is a significant increase in average annual discharge of 57.96 mm (p = 1e-5; one-sample t-test) relative to the baseline period (Figure 5).Of the mean overall change, climate is a slightly stronger contributor than land use, though both have significant impacts.Climate change causes a mean 36.94 mm increase in discharge (p=0.002;63.7% of overall change), while land use change contributes a 21.02 mm increase (p=0.002;36.3% of overall change) relative to the baseline period.However, there is substantial interannual variability in the relative strength of the two drivers.Overall changes in discharge relative to the baseline period appears to respond most strongly to climate variability from year-to-year, with a consistent but low-level positive effect due to land use change (Figure 5d).Both climate and land use effects are positive in 17 of 21 years (Figure 5e).
Quickflow and baseflow contribute approximately equally to the observed increases in discharge, with a mean annual increase in baseflow of 28.18 mm (p<0.0001; Figure S5) and mean annual increase in quickflow of 29.79 mm yr -1 (p=3e-5; Figure S6).However, the relative contribution of land use and climate to these two components of overall discharge varies.For quickflow, the increase is dominated by climate (23.53 mm; p=0.002) with a small contribution from land use (6.27 mm; p=0.06).In contrast, for baseflow the increase due to land use is larger (15.51 mm; p=6e-5) than the increase due to climate (11.26 mm; p=0.007); however, the proportion of total change in baseflow attributed to land use may be an overestimate due to long timescales of baseflow response to changes in watershed-scale subsurface storage (see discussion in Section 4.1).
Using our continuous PLSR-based approach, we also identify changes through time in the relative contribution of climate and land use during the prediction period.While discharge is increasing through time at a rate of 1.2 mm yr -1 , this trend is not significant (p=0.2).This trend is dominated by land use effects (0.79 mm yr -1 ; p=0.13), while climate effects are smaller (0.38 mm yr -1 ; p=0.73).The prediction period corresponds with a significant positive trend in the percent of the watershed with urban land use (1.18 %/year; p<0.05).The trend in discharge is primarily driven by increases in baseflow, which has positive overall (1.80 mm/yr; p=0.02) and land use trends (1.15 mm yr -1 ; p=0.02) with no significant climate trend (p=0.31)during the prediction period (Figure S5).In contrast, there are no significant quickflow trends for overall, land use, or climate changes (Figure S6).

Sensitivity analysis of baseline period
While the results described above all use a baseline period of 1974-1995, model performance is good (NSE>0.59)and comparable regardless of the baseline period used (Figure 6a).Validation performance does increase notably (NSE>0.68)when 1993, a particularly wet year, is included in the baseline period, and the highest model skill occurs with a baseline period end year of 1997 (NSE=0.761).Similarly, the relative importance of land use and climate are comparable for all baseline periods at both a mean and interannual scale (Figure 6b-c).Comparing within the common prediction period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), there are no significant differences between baseline period end years in PLSR-estimated changes due to either land use or climate change (Figure 6b-c).

Sensitivity analysis of input meteorological data
Our results also indicate that the attribution of changes to land use vs. climate change is relatively insensitive to the choice of input meteorological data (Figure 7).The PLSR models built using either meteorological input dataset are able to accurately predict discharge during the baseline period, though performance is better with Airport meteorological data (NSE=0.839)than Arlington meteorological data (NSE=0.695).Predictions from both input datasets indicate the same direction of mean annual changes in discharge due to climate for 19 of the 21 years during the prediction period, with a strong positive correlation (r = 0.91) between the two methods.Furthermore, a two-sample Kolmogorov-Smirnov test indicates that there is not a statistically significant difference between the distributions of changes due to climate during the prediction period from the two data sources (p=0.987).Overall, differences between PLSR results from the two input meteorological sources were small, with an annual root mean squared difference of 25.1 mm (13.4% of total interannual variability).
3.2 Future scenario analysis

Model validation
When analyzing output from the simulated future scenarios, monthly PLSR models perform very well, with NSE of 0.983, 0.870, and 0.963 for ET, drainage, and direct runoff, respectively (Figure 8).RMSE are 4.60 mm (3.6% of range of observations), 4.09 mm (5.3%), and 2.30 mm (2.8 %), respectively.Performance is also strong at an annual level, with NSE of 0.843, 0.958, and 0.967 for ET, drainage, and direct runoff.Statistics summarizing overall and monthly fits for each hydrological flux are provided in Table S1.

Climate and land use impacts on the water balance
The scenarios generated a wide range of climate and land use model inputs that exposed the relative impacts of these two drivers under a variety of conditions, with AI and AR representing the extremes for most inputs (Figure 3).For example, while air temperature increased in all scenarios relative to the historical period, there is ~4°C difference across the four scenarios, with the most extreme increase in AR and the mildest increase in AI.Similarly, precipitation changes varied across the four scenarios, with ~400 mm of variability between scenarios; AR had the most extreme increases in precipitation, particularly during the 2030s.Land use change also varied substantially between scenarios; row-crop agriculture, for example, was relatively consistent through time in the AI scenario, but decreased in each of the other scenarios and was almost completely eliminated by the end of the AR scenario.Urban land use was highest for the AI scenario, lowest in the AR scenario, and relatively unaffected in the CC and NW scenarios.
Overall changes in watershed-average ET are uniformly positive relative to the baseline period across all combinations of scenarios, ranging from 23.49 mm yr -1 to 90.83 mm yr -1 over the final two decades of the simulations (Figures 9a, 10a).These increases are dominated by climate effects (24.77 mm yr -1 to 70.71 mm yr -1 ), with a small but typically positive effect of land use (-1.55 mm yr -1 to +20.12 mm yr -1 ).While climate and land use have synergistic effects across most scenarios, the effects of land use most strongly counteract those of climate in the AR scenario.AR is characterized by decreases in row-crop agriculture and increases in natural vegetation, while land use effects are closest to 0 in the AI scenario, which is characterized by widespread expansion of impervious cover.Patterns in ET through time correspond primarily to changes in temperature and reference ET.For example, in all scenarios with AI climate ET peaks in the 2040s, declines through the 2050s to a low in ~2060, and rises in the final decade of the simulations (Figure 9a); this pattern corresponds with temperature in the AI scenario, which is one of the primary controls on reference ET (Figure 3).Variability in the predicted land use and climate effects increase with time through the scenario as climate deviates further from the baseline climate used to develop the PLSR relationships, indicating that our permutation-based approach is able to adequately identify and express uncertainty.
There is more temporal variability in drainage results compared to ET, with overall mean changes ranging from -142.12 mm yr -1 to 65.17 mm yr -1 over the final two decades of the simulations (Figures 9b, 10b).Both climate and land use can have positive effects (increase in drainage) and negative effects (decrease in drainage), with interannual variability in drainage corresponding closely to climate.Climate effects range from -74.58 mm yr -1 to 38.00 mm yr -1 , and land use effects from -67.55 mm yr -1 to 27.17 mm yr -1 .Unlike ET, however, climate-driven and land use-driven changes do not have a consistent synergistic or antagonistic relationship, with a primarily synergistic interaction in the AI and AR climate scenarios and a primarily antagonistic interaction in the CC and NW climate scenarios (Figure 10b).However, this directional change is not constant through time.Across all scenarios with AR land use, in particular with AI and AR climate, the effects of land use on drainage are fairly neutral through the 2030s then positive in the 2040s and 2050s (Figure 9b), a period characterized by a decrease in urban land use and increase in natural vegetation (Figure 3).While ET seems to be driven primarily by temperature, changes in drainage respond more to the relative balance of ET and precipitation.In the AR climate scenario, changes in drainage relative to the baseline period begin declining from their peak in the late 2040s, becoming negative in the mid-2050s and plateauing in around 2060 for the remainder of the simulation.This trend coincides with a period of decreasing precipitation and increasing reference ET (Figure 3).
Like ET, direct runoff increases in all future scenarios, with increases ranging from 8.05 mm yr -1 to 49.20 mm yr -1 over the final two decades of the scenarios (Figures 9c, 10c).As with ET, the effects of climate tend to dominate with land use effects mostly contributing to a small but synergistic effect: climate accounts for -1.24 mm yr -1 to +39.79 mm yr -1 of overall changes, compared to 7.20 mm yr -1 to 11.47 mm yr -1 for land use.Through time, changes in direct runoff track total annual precipitation, annual extreme precipitation events, and annual reference ET (Figure 3a,b,d).For example, in the AR climate scenarios, changes in direct runoff are largest in the 2030s and 2040s (Figure 9c), the wettest period on record which included the largest number of extreme precipitation days (Figure 3).In contrast, the NW climate scenarios have a decline in direct runoff from the 2040s through the end of the simulation (Figure 9c) which occurs despite increasing overall and extreme precipitation due to increasing temperature and reference ET (Figure 3).

Historical changes in discharge
Our results for the PBS indicate that land use contributes to ~40% of observed increases in discharge while climate contributes ~60%, shedding light on a trend of unknown origin previously documented by Gebert et al. (2012).Furthermore, our method's ability to provide continuous results through time reveals the changing relative importance of these two drivers: land use-induced changes in discharge are increasing through time at approximately twice the rate of climate-driven changes.This trend appears to be driven primarily by a strong trend of increasing urban land cover, with a land use-driven increase in discharge of 2.12 mm for each 1% increase in urban land use within the PBS.
Disentangling these two drivers, as well as their changes through time, provides insight into the effects of historical watershed-scale management decisions.While urbanization-driven increases in discharge are often associated with increased runoff (Boggs & Sun, 2011;Debbage & Shepherd, 2018;Rose & Peters, 2001), our results indicate that increases in quickflow and baseflow are comparable at the monthly scale (our analysis is done at a monthly timestep and is not intended to capture effects at the event scale).Moreover, overall and land use effects on baseflow are increasing through time, unlike quickflow.
Given that the period of study coincides with an expansion of urban and impervious cover, the significant effect of land use change on baseflow, not quickflow, is surprising.While outside the scope of the present study, we suggest two possible explanations which may be operating in tandem.First, the observed increase in baseflow may demonstrate that strict infiltration requirements for new developments in the PBS (Ch.26.06(3),City of Middleton ordinances) are successfully reducing the impacts of climate change and urbanization on direct runoff, but are increasing groundwater recharge and baseflow due to more focused infiltration as well as other potential water sources associated with urbanization (e.g.urban irrigation and leaky water infrastructure).Second, our PLS-based methodology may be attributing the effects of long-term increases in groundwater storage to land use change.There is a long-term increasing trend in groundwater levels of 0.3 m decade -1 with substantial variability at yearly to decadal timescales and a nonlinear response of baseflow to water table depth (Figure S7).Since changes in storage are endogenous to the PBS and the timescale over which groundwater storage changes are longer than the maximum timescale considered in our PLSR relationships (one year), baseflow response to changes in watershed-scale storage could be methodologically attributed to the effects of land use change, which tends to follow long-term trends but has little interannual variability.
Combined, these results may help guide future management interventions targeted at buffering the observed changes in discharge, quickflow, and baseflow associated with urbanization.Infiltration-based stormwater management (e.g.distributed green infrastructure) may have an unintended effect of increasing baseflow, potentially creating more drought-resistant streams.Given that infiltration-based stormwater management is also effective at counteracting climate-induced changes in discharge during extreme events, these practices may present an opportunity to protect aquatic ecosystems during both low-and high-flow periods.However, research elsewhere has found that reductions in runoff volumes do not always translate to increased baseflow due to watershed-specific factors such as the amount and distribution of impervious cover (Fanelli et al., 2017) and enhancing groundwater recharge using green infrastructure may lower water quality in surficial aquifers (Andres et al., 2018).This highlights a need to better understand how land use propagates through groundwater flow systems to impact downstream terrestrial and aquatic ecosystems (Bhaskar et al., 2016(Bhaskar et al., , 2018;;Breyer et al., 2018;Jefferson et al., 2017;Zipper et al., 2017a).

Future scenario analysis
Results from our factorial set of scenarios indicate that the effects of climate, not land use change, will likely dominate the future water balance of the YW.Specifically, ET seems to respond most strongly to temperature, while direct runoff responds most strongly to precipitation.Climate effects on drainage are driven primarily by the balance of supply (precipitation) and demand (reference ET).As precipitation projections have considerably less certainty than temperature projections (WICCI, 2011), this makes understanding the impacts of climate and land use change on surface water and groundwater resources particularly challenging.In fact, the similarity of predicted land use effects between different land use scenarios for a given climate (e.g.columns in Figures 9 and 10) indicates that errors in the PLSR relationships due to climate change may be larger than the effects of land use change.
While the effects of land use are smaller than those of climate, several key patterns and interactions with climate emerge.Fluxes occurring at the land surface (ET and direct runoff) tend to have synergistic relationships between climate and land use change, with increases resulting from climate change constituting the majority of overall change.In contrast, drainage has a mix of synergistic and antagonistic relationships and the largest land use effects of any of the fluxes studied, exceeding 50% of overall change in some combinations of land use and climate scenarios.This indicates that, while the effects are relatively small, land use changes can act as a buffer on drainage from climate change at a watershed scale.While groundwater recharge is typically thought of as beneficial, excess groundwater can have negative effects on several ecosystem services including reductions in flood retention capabilities, risk of basement flooding in urban areas, and decreases in agricultural productivity associated with oxygen stress (Booth et al., 2016a).It is therefore critical to consider the implications of either an increase or decrease in watershed-scale drainage for groundwater flow and associated ecosystems when making land use decisions.
Additionally, the AI land use scenario (characterized by urbanization) and AR land use scenario (characterized by a return to natural ecosystems) consistently have the most extreme impacts on the water balance.Across all scenarios, AI has the smallest effect on ET and drainage, but the largest effect (most negative) on direct runoff.In contrast, AR has the largest effect on drainage (most positive) and ET (most negative), and among the smallest effects on direct runoff.This highlights the important role of land use in determining the partitioning of water at the land surface and in the root zone.

Synthesis and management implications
Both the historical discharge analysis in the PBS and the future scenario analysis of the YW indicate that climate is the key control over the water balance, though the analyses differ in the relative importance of land use.In the PBS, results indicate that climate change contributes to ~60% of observed changes in discharge, with approximately equal impacts on quickflow and baseflow, though the effects of land use are increasing through time.Similarly, the simulated future scenario analysis points to climate as the key control over direct runoff (as well as ET and drainage), with relatively smaller effects of land use.To better assess potential causes of these differences, we extracted Agro-IBIS direct runoff output from the Pheasant Branch portion of the YW and repeated all analyses for the common period of record .While the baseline period data differs between the two analyses due to different meteorological input data in the 1974-1985 period (see section 2.2.3.2),results for the overall degree of change are comparable.Results from historical analysis of the Agro-IBIS output for Pheasant Branch (Figure S8) finds that 52.3% of the overall changes in direct runoff during the prediction period result from climate and 47.7% result from land use (compared to 79.0% climate and 21.0% land use for quickflow estimated from baseflow separation; Figure S6).Also similar to the results from baseflow separation, there is no significant trend through time for overall, land use, or climate-induced changes in Agro-IBIS direct runoff for the portion of the prediction period with real climate inputs (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013).This analysis indicates that the differences between the historical discharge analysis and the future scenario analysis is driven by several factors.First, the degree and trajectory of land use change varies between the spatial scales used for the two analyses.The PBS is significantly smaller than the YW (~3% of the YW) and has experienced relatively unidirectional land use change (urbanization) during the historical period (Figure 5c).In contrast, the future scenarios include a large variety of contrasting land use changes which may partially counteract each other when aggregated to the watershed scale.The stronger land use signal in the PBS relative to the YW implies that, just as the impacts of climate change on streamflow are attenuated in larger river networks (Chezik et al., 2017), so too can larger spatial scales attenuate the effects of land use change.Second, climate change during the future scenario analysis (2°C to 5.5°C warming) is more extreme than has been observed in the historical record.Third, our biophysical modelling approach may underestimate differences in hydrological properties between land uses (see Section 4.4).
While our analysis agrees with recent work showing that climate effects may dominate future hydrological changes (Martin et al., 2017;Peng et al., 2016;Pribulick et al., 2016;Wang et al., 2017), we also highlight the critical need to target land use interventions locally to maximize benefits in areas of concern (Fry & Maxwell, 2017;Schifman et al., 2017).The results presented for the YW average hydrologic response over an area of 1344 km 2 , and therefore neglect spatial heterogeneity in land use which can impact the local water cycle (Deshmukh & Singh, 2016;Frans et al., 2013;Haddeland et al., 2007;Wang & Stephenson, 2018).As observed in the historical discharge analysis, management interventions can impact hydrological processes at a nearly comparable level to climate change.
Elsewhere, previous work has shown that, for example, the expansion of biofuel cropping systems can change ET (Harding et al., 2016;Joo et al., 2017;VanLoocke et al., 2010;Wagner et al., 2017); land use change can either reduce or increase groundwater recharge (Giménez et al., 2016;Newcomer et al., 2014;Oliveira et al., 2017;Qiu & Turner, 2015;Robertson et al., 2017;Zipper et al., 2017a); and urban green infrastructure and agricultural drainage management can successfully reduce runoff volumes (Allred et al., 2003;Elliott et al., 2016;Schott et al., 2017;Shuster et al., 2017;Wadzuk et al., 2010).Each of these represents a management practice that can alter a hydrologic flux of interest in the context of climate change which may have significant local benefits.

Methodological strengths and limitations
While statistical regression techniques have previously been used to separate the impacts of climate and land use on streamflow (Section 2.1), our technique has several novel contributions.First, users of regression-based methods typically divide their data into two discrete chunks ("before" and "after") and separate the temporally-averaged land use and climate impacts using residuals from the "after" period.In reality, of course, both land use and climate change are rarely step changes, but rather shift gradually over time.While the approach introduced here uses a baseline period, it also provides continuous estimates of the relative importance of land use and climate change over time during both the baseline and prediction period, which makes it possible to identify trends in drivers of hydrological change.For example, in the PBS, we reveal that the impacts of land use change are increasing over time for both discharge and baseflow; while climate change has significantly increased streamflow but there is no significant trend during the prediction period.Second, the continuous separation through time makes it possible to assess synergistic and antagonistic relationships between climate and land use, as well as the changing nature of these interactions through time.Third, as opposed to MLR which is frequently used to separate land use and climate change effects (Huo et al., 2008;Jiang et al., 2011;Ye et al., 2003;Zhang et al., 2016), we use a PLSR model and evaluate its performance relative to MLR and PCR model.Our PLSR approach relies on automated significance-pruning to select predictor variables from a set of candidates, thus reducing potential spurious correlations and potential researcher biases, and eliminated issues associated with collinearity to provide more robust predictions (Hwang & Nettleton, 2003).
We do, however, note several potential limitations to our method, most of which are common to all statistical approaches for separating climate and land use change.For instance, developing statistical relationships based on one period of time and applying them to another may result in extrapolation beyond the conditions for which the relationships are well-suited.This problem is common to all regression-based methodologies and may be particularly challenging in the context of nonstationarity or where emergent properties of the relationship between climate and land use change lead to novel future responses (Milly et al., 2008), or where significant changes in watershed storage occur (e.g.rising groundwater levels as discussed in Section 4.1).In our case, sensitivity analysis results demonstrate that separation of climate and land use effects is relatively insensitive to the selection of the baseline period and meteorological input data, as long as the performance of the PLSR model is validated and demonstrated to accurately reproduce observations, thus minimizing concerns regarding nonstationarity.However, alternate approaches to selecting model calibration periods may provide a useful approach to baseline period selection (Razavi & Tolson, 2013).Additionally, using a permutation-based approach to fitting models allows us to identify and visualize increases in uncertainty during climate conditions outside the range of the baseline period (e.g. Figure 9).Additionally, regression-based approaches assume that climate and land use have linearly additive effects (e.g.Eq. 4); however, due to the complexity of socio-environmental systems, nonlinear hydrological responses to climate and land use change may violate this assumption (Krause et al., 2017).Thus, regression-based methodologies may be less accurate when nonlinear behaviors are present, and process-based models may be more effective (Pribulick et al., 2016;Pumo et al., 2017).Furthermore, while Agro-IBIS is a state-of-the-art dynamic vegetation and agroecosystem model, some land use characteristics which may impact the water cycle are not represented.For example, changes in soil hydraulic properties between land uses and through time are not simulated (Paturel et al., 2017), nor are soil hydraulic properties coupled to soil organic content (Ankenbauer & Loheide, 2017).Improving parameterizations and including these processes would likely increase the simulated differences between land use types in the future scenario analysis and increase the relative importance of land use.Our statistical relationships also do not take into account other factors which may drive changes in the water balance; for example, each scenario has a representative atmospheric CO2 concentration pathway (Booth et al., 2016b).Given that carbon and water cycles are coupled in Agro-IBIS via stomatal conductance, CO2 may also be a relevant predictor variable, particularly for ET and under conditions with significant land use change between C3 to C4 vegetation (Twine et al., 2013).While our methodological framework is sufficiently flexible to include variables such as CO2, we elected to exclude it from analysis in order demonstrate our methodology using easily obtained meteorological data.

Conclusions
This study introduces a new methodological framework to separate the effects of climate and land use on the water cycle continuously through time.We tested this approach using three different regression models (PLSR, PCR, and MLR) and found the best performance with PLSR.We then applied our approach to both observed and modeled data for the YW in south-central Wisconsin.Analysis of historical discharge data for the PBS indicated that climate change has caused ~60% of the observed changes in discharge over the past two decades, with a significantly increasing impact of land use change (urbanization) on both baseflow and overall discharge.Using a factorial combination of four contrasting land use and climate scenarios, we show that future changes in the YW's land surface water balance (ET, drainage, and direct runoff) are likely to be dominated by effects of climate change: ET was most affected by changes in temperature, direct runoff by changes in precipitation, and drainage by changes in both precipitation and reference ET.Land use effects were larger on drainage than either ET or direct runoff.
Overall, these results indicate that the effects of land use and climate are not static through time, and separating the relative contribution of these two variables to hydrological change should not be done via the simple separation of time into discrete periods; rather, it must be done in a continuous manner.Furthermore, we show that using land use to mitigate the effects of climate change on the water cycle may be challenging in large watersheds which contain a diversity of land use trajectories.However, our results indicate that the effects of land use change are larger in the PBS than the YW as a whole due to the relatively monodirectional land use change from agriculture to urbanization.Therefore, local management interventions targeted at subwatershed scales to achieve specific desired outcomes may be an effective path forward to protecting water resources from future climate change.(Bristow & Campbell, 1984).(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)) (Fry et al., 2011;Homer et al., 2007Homer et al., , 2015));

Figures
Figures Figure 1.Flowchart illustrating method for separating land use and climate effects, demonstrating the monthly resolution used in the current study.

Figure 2 .
Figure 2. (a) Map of Yahara Watershed showing land use in 2014 at Agro-IBIS model resolution (220 m grid cells).The green dot shows the Pheasant Branch gauging station, and thick black line outlines the contributing area.(b) Relative proportion of different land uses in the Yahara River Watershed and Pheasant Branch Subwatershed.The Agriculture class includes all crops and pasture (top 7 legend entries in panel a).The Urban class includes all urban density levels as well as barren land.Natural includes forest, grassland, and wetlands.

Figure 3 .
Figure 3. Annual watershed-average meteorological (a-d) and land use (e-h) model input for the four scenarios.In each plot, the gray shading represents the historical (1986-2013) range.Meteorological variables (a-d) are smoothed with an 11-year moving average.Plot show (a) annual cumulative precipitation; (b) annual days with >1" (25.4 mm) precipitation; (c) mean maximum daily temperature; (d) mean Penman-Monteith reference evapotranspiration; (e) percent of domain with corn land cover; (f) percent of domain with urban (low, medium, and high density) land cover; (g) percent of domain with deciduous forest land cover; (h) percent of domain with wetland land cover.

Figure 4 .
Figure 4. Comparison of difference statistical models: (a-b) partial least squares regression (PLSR); (c-d) principal components regression (PCR); and (e-f) multiple linear regression (MLR) for the Pheasant Branch Subwatershed.Left column shows monthly mean of all validation samples (colored line) compared to measured discharge (black line).Right column of plots shows annual measured discharge (black line) and mean +/-1 standard deviation of predicted discharge (color ribbon).(g) shows monthly distributions of observations (white box) and validation samples (colored boxes) for each method.

Figure 5 .
Figure 5. Results from analysis of Pheasant Branch historical discharge data.(a) Comparison between observed and predicted (mean of random validation samples for all PLSR permutations) for baseline period; (b) boxplots showing monthly distributions of discharge for observed (all years) and predicted (all years and all permutations); (c) percent of Pheasant Branch Watershed with urban land use (combined high, medium, and low density) from WISCLAND (Wisconsin Department of Natural Resources, 2016) and NLCD datasets(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)); (d) change relative to baseline period, with solid line showing overall change and ribbons spanning +/-1 standard deviation of the mean across all permutations; (e) density plot of mean annual changes in discharge due to land use, climate, and overall.Legend in (a) also applies to (b) and legend in (d) also applies to (e).

Figure 6 .
Figure 6.Sensitivity of PLSR results for Pheasant Branch watershed to selection of baseline period.(a) Nash-Sutcliffe Efficiency for the calibration and validation samples as a function of the end of the baseline period.Changes in discharge due to (b) land use and (c) climate, color-coded by the end of the baseline period.All baseline periods begin in 1974.For results in Figure 5 and discussed in text, baseline period ends in 1995 (black line).

Figure 7 .
Figure 7. Sensitivity of attributed changes in Pheasant Branch discharge due to (a) land use change and (b) climate change to choice of meteorological input data for generating PLSR relationships.Ribbons show mean +/-1 standard deviation from all permutations.

Figure 10 .
Figure 10.Density plots showing distribution of overall (black line), climate-induced (shaded red), and land use-induced (shaded green) changes to annual watershed-average (a) evapotranspiration, (b) drainage, and (c) direct runoff.Distributions are for the final 20 years of the simulation (2051-2070), relative to the baseline period (1986-2013).In each set of 16 plots, the labels along the top show the climate scenario and the labels along the right side show the land use scenario.

Figure S2 .
Figure S2.Comparison between (a) daily mean and (b) monthly mean incoming shortwave radiation.Arlington data are from agricultural weather station at Arlington Agricultural Research Station.GHCN data are estimated using daily maximum and minimum temperatures from a calibrated Bristow-Campbell equation(Bristow & Campbell, 1984).
Figure S5 (as Figure 5, but for baseflow).Results from analysis of Pheasant Branch historical baseflow data.(a) Comparison between observed and predicted (mean of random validation samples for all PLSR permutations) for baseline period; (b) boxplots showing monthly distributions of baseflow for observed (all years) and predicted (all years and all permutations); (c) percent of Pheasant Branch Watershed with urban land use (combined high, medium, and low density) from WISCLAND (Wisconsin Department of Natural Resources, 2016) and NLCD datasets(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)); (d) change relative to baseline period, with solid line showing overall change and ribbons spanning +/-1 standard deviation of the mean across all permutations; (e) density plot of mean annual changes in baseflow due to land use, climate, and overall.Legend in (a) also applies to (b) and legend in (d) also applies to (e).
Figure S5 (as Figure 5, but for baseflow).Results from analysis of Pheasant Branch historical baseflow data.(a) Comparison between observed and predicted (mean of random validation samples for all PLSR permutations) for baseline period; (b) boxplots showing monthly distributions of baseflow for observed (all years) and predicted (all years and all permutations); (c) percent of Pheasant Branch Watershed with urban land use (combined high, medium, and low density) from WISCLAND (Wisconsin Department of Natural Resources, 2016) and NLCD datasets(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)); (d) change relative to baseline period, with solid line showing overall change and ribbons spanning +/-1 standard deviation of the mean across all permutations; (e) density plot of mean annual changes in baseflow due to land use, climate, and overall.Legend in (a) also applies to (b) and legend in (d) also applies to (e).
Figure S6 (as Figure 5, but for quickflow).Results from analysis of Pheasant Branch historical quickflow data.(a) Comparison between observed and predicted (mean of random validation samples for all PLSR permutations) for baseline period; (b) boxplots showing monthly distributions of quickflow for observed (all years) and predicted (all years and all permutations); (c) percent of Pheasant Branch Watershed with urban land use (combined high, medium, and low density) from WISCLAND (Wisconsin Department of Natural Resources, 2016) and NLCD datasets(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)); (d) change relative to baseline period, with solid line showing overall change and ribbons spanning +/-1 standard deviation of the mean across all permutations; (e) density plot of mean annual changes in quickflow due to land use, climate, and overall.Legend in (a) also applies to (b) and legend in (d) also applies to (e).
Figure S6 (as Figure 5, but for quickflow).Results from analysis of Pheasant Branch historical quickflow data.(a) Comparison between observed and predicted (mean of random validation samples for all PLSR permutations) for baseline period; (b) boxplots showing monthly distributions of quickflow for observed (all years) and predicted (all years and all permutations); (c) percent of Pheasant Branch Watershed with urban land use (combined high, medium, and low density) from WISCLAND (Wisconsin Department of Natural Resources, 2016) and NLCD datasets(Fry et al., 2011;Homer et al., 2007Homer et al., , 2015)); (d) change relative to baseline period, with solid line showing overall change and ribbons spanning +/-1 standard deviation of the mean across all permutations; (e) density plot of mean annual changes in quickflow due to land use, climate, and overall.Legend in (a) also applies to (b) and legend in (d) also applies to (e).
Figure S7.(a) Water table depth at a well just north of the Pheasant Branch subwatershed in an area that has seen minimal urbanization (USGS 430638089353101); (b) estimated monthly baseflow in the Pheasant Branch watershed; and (c) baseflow as a function of water table depth.Blue lines in (a-b) show linear best-fit with 95% confidence interval.Red line in (c) shows the lower bound of baseflow increasing nonlinearly as a function of water table depth.

Figure S9 .
Figure S9.Comparison of simulated and observed direct runoff at (a) monthly and (b) annual timescales for the Pheasant Branch subwatershed.Observed runoff is quickflow estimated from overall discharge at the Pheasant Branch gauging station using baseflow separation.Agro-IBIS runoff is the average of direct runoff at all grid cells within the contributing area of the Pheasant Branch gauging station.