Relative impacts of climate change and land cover change on streamflow using SWAT in the Clackamas River Watershed, USA

Global climate change (CC) and rapid urbanization are likely to have strong impacts on water resources around the world (IPCC 2014). Water scarcity, distribution, and access to water remain as some of the biggest challenges in the 21st century. Billions of people globally will not have sustainable access to clean drinking water due to the impacts of global warming (Mukheibir 2010; Schewe et al. 2014). CC impacts the hydrologic cycles across multiple scales (Arnell & Gosling 2013; Hattermann et al. 2017). Watershed scale hydrological predictions rely on the transfer of large-scale climate variables to more regional meteorological factors such as precipitation and temperature. The multi-ensemble means of different general circulation models (GCMs) have been popular among researchers for projecting future climate and impacts on streamflow. However, using multiple GCMs as inputs may increase data and modeling uncertainties, as climate and water resource projections vary between each GCM (Guimberteau et al. 2017; Thompson et al. 2017; Shen et al. 2018).

Under CC, the Willamette River Basin (WRB) in the Pacific Northwest (PNW) region will exhibit significant changes in water balance and temperature (Jaeger et al. 2017). The magnitude of change will vary based on seasonality and location, as well as regional climate interactions with land cover (LC) and land use (Jung & Chang 2011; Vano et al. 2015b). Catchments in the WRB rely heavily on snowpack for summer water supply. The projected change in precipitation patterns showed that more precipitation would fall as rain than snow and snowpack will be drastically reduced in the Cascade range (Catalano et al. 2019). Snow-fed rivers in the WRB provide essential water resources for irrigation and municipal consumption. Under CC, short-term drought risk is projected to increase in the summer due to earlier snowmelt and less precipitation (Jung & Chang 2012). As population increases, urban development sprawls toward city boundaries and converts rural landscapes into more impervious landscapes (Hoyer & Chang 2014). Water demand grows with population increases, adding stress to the currently vulnerable water system that is impacted by recent extreme climatic events in the region such as the 2015 drought (Marlier et al. 2017). Water demand in the Portland metropolitan area of the WRB is expected to increase in the coming decades under GCM scenarios and projected land-use change (Parandvash & Chang 2016).

Hydrologic modeling using data from downscales GCMs have underlying uncertainties that are yet to be quantified and resolved (Jung et al. 2012; Hattermann et al. 2018; Her et al. 2019). As trends of CC and urban development continue throughout the 21st century, researchers need to improve modeling techniques to more accurately predict the combined effects of climate and land use on water quantity and quality (Praskievicz & Chang 2009; Xie & Lian 2013; Chen et al. 2019). Models such as the Better Assessment Science Integrating Point and Nonpoint Sources (BASINS) and Windows-based Hydrologic Simulation Program (WinHSPF) have also been used to access the separate and combined impacts of climate and land cover (CLC) change scenarios in surrounding watersheds (Praskievicz & Chang 2011). The U.S. Geological Survey's (USGS) Precipitation–Runoff Modeling System model was run by researchers in the region to model changes in runoff, and hydrological variability was expected to increase with seasonal flow becoming more sensitive to CC (Chang & Jung 2010). The close relationship between CC and land cover change (LCC) calls for a more systemic approach to modeling hydrology (Sterling et al. 2013; Devia et al. 2015; Dwarakish & Ganasri 2015). The semi-distributed Soil and Water Assessment Tool (SWAT) model has been a useful tool utilized frequently by researchers to predict runoff, nutrient, and sediment transport (Raneesh & Santosh 2011; Arnold et al. 2012; Leta et al. 2016; Suttles et al. 2018; Hajihosseini et al. 2019). Additionally, SWAT is well capable of simulating and projecting hydrologic responses to CC and LCC in small to large watersheds, by allowing the feedback between CC and evaporative demand (Praskievicz & Chang 2009; Kim et al. 2013; Psaris 2014; Tan et al. 2014; Pervez & Henebry 2015).

This study aims to investigate the hydrologic response to CC and LCC of the Clackamas River Watershed (CRW) in the WRB. CC and LCC closely interact with each other on multiple scales when used to predict hydrologic responses (Lahmer et al. 2001). The separate and combined impacts of CLC change on hydrology (Zhang et al. 2018) have been researched and applied in projecting streamflow (Kim et al. 2013; Zhang et al. 2015), stormwater runoff (Talib & Randhir 2017), water availability (López-Moreno et al. 2014), and water quality (Sun et al. 2016). However, the development of LCC scenarios that are representative of the variability and range in climate modeling remains difficult. Past studies have used statistically downscaled GCMs to model LCC, but uncertainties can still occur during data processing and scaling (Solecki & Oliveri 2004; Tan et al. 2014; Prestele et al. 2016). To reduce assumptions and uncertainties in hydrologic models, LCC modeling efforts must include a subset of CC models that are representative of a range of possible future scenarios (Turner et al. 2015; Vano et al. 2015a).

In terms of scale, past studies (He et al. 2013; Farinosi et al. 2019) often modeled and compared hydrologic responses at the basin scale, between sub-basins, or across two different basins with similar climate and topography. As humans influence CLCs with increasing urbanization, modeling CLC changes on a landscape gradient from urban to rural is becoming more critical, as it can distinguish land-use change processes and different types of disturbances across the landscape (Clavero et al. 2011). The goal of this study is to understand the separate and combined impacts of CLC change through streamflow indices that are able to represent the timing and magnitude of flow over time. Low flow, peak flow, seasonal mean flow, and center timing (CT) of flow are all useful streamflow indices used to predict spatial and temporal changes in runoff of complex watershed systems facing CC and urbanization (Chang & Jung 2010; Choi et al. 2017). Unlike previous studies, our work aims to use these streamflow metrics as indicators of change and also incorporate tightly coupled CLC change scenarios (individual and combined) into our models to yield a reasonable range of impact scenarios between the lower and upper watershed in the near future (2050s) and the distant future (2080s).

We chose the CRW, located geographically between longitudes 121° 45′ 12″ and 122° 36′ 25″ E and latitudes 44° 49′ 26″ and 45° 22′ 20″ N in the Lower WRB of Oregon in the United States as the study area (Figure 1). The lower part of the watershed is heavily urbanized with medium- to low-density developments. Geology of the watershed is dominated by the western cascade volcanic rocks, with a small portion of the lower watershed falling in the Willamette Valley alluvium deposits and the very upper part of the watershed falling inside the high cascade range with colder and more snow deposits. The total population of the watershed is approximately 200,000 people across a drainage area of 2,435 km². The entire watershed consists of 5% urban developments, 10% agricultural lands, and 85% forested lands. However, urban developments and agricultural lands are more concentrated on the lower watershed, while forested lands dominate the middle and upper parts of the watershed. The main stem of the Clackamas River originates from Mount Hood, flowing through the pristine mixed forest from southeast to northwest into the Willamette River. The Clackamas River flows through both rural and urban areas as well, providing drinking water to roughly 350,000 people within and adjacent to the watershed. Drinking and wastewater treatment plants are all located near the mouth of the river downstream as well as a USGS stream gauge with continuous discharge and water quality monitoring. The study area is highly vulnerable to CC, as it is heavily dependent on diminishing snowpack for water supply and highly sensitive to wet season floods and dry season droughts (Graves & Chang 2007). Both high and low flows are concerns for water managers, as high flows typically accompany turbid water (Chen & Chang 2019), while low flows reduce available water for drinking and irrigation.

The climate of the study area is considered Mediterranean, with a prolonged winter and fall rainfall period and dry, warm summer. Climate data from 1981 to 2010 showed that mean air temperature is approximately 4.6 °C in January and 20 °C in July. Precipitation is the most abundant during December, averaging 183 mm and driest in July with only 19 mm. The majority of precipitation falls as snow on the upper part of the watershed and becomes essential runoff in the following spring and summer. The mean annual runoff is 138 cm/year in the watershed from 1981 to 2010. Runoff patterns vary by season with highs during late winter early spring and lows in mid-summer. With a growing population and drinking water demand, LC in the watershed is expecting significant changes in the year 2040–2070 and 2070–2100 based on projected CC and LC (Turner et al. 2015). We divided the watershed into the lower and upper watersheds. The lower watershed represents the more urbanized and agricultural heavy area, while the upper watershed consists of mainly evergreen forest.

SWAT was selected to model hydrologic changes under individual and combined CLC scenarios in our studied watershed due to its ability to capture physical hydrology processes at the watershed to basin-level on a continuous-time scale. The SWAT is a semi-distributed continuous-time model capable of modeling streamflow, sediment transport, and nutrient runoff at the watershed to basin scale (Arnold et al. 2012). Through 30 years of research and development, SWAT is well-documented with a user-friendly interface. Unlike other process-based hydrologic models, SWAT has its own calibration and sensitivity analysis software called SWAT-CUP and is capable of running on geospatial platforms like ArcGIS and QGIS. Temporally SWAT is capable of simulating past to future periods from the scale of sub-hourly to a century (Douglas-Mankin et al. 2010). Spatially being semi-distributed, SWAT can model hydrologic cycle using physical-based water balance equations at three different scales from the entire watershed, subwatershed, to lastly hydrologic response units (HRUs). HRUs are unique pixels below the subwatershed level that are combinations of LC, soil type, and slope derived from digital elevation models (DEMs). This research was done using the ArcSWAT plug-in in ArcMap 10.5.1 (ESRI 2017) where DEM input was also essential in delineating watershed boundaries and outlets. LC data and climate data are also required as inputs to run the model on the desired period. Hydrologic responses are simulated at the HRU level then aggregated to show streamflow at the subwatershed and watershed scale (Arnold et al. 1998). Details of all the required input data sources and description are listed in Table 1. Evapotranspiration (ET) is an important component in the water balance equation when estimating streamflow, and we used the Hargreaves method for estimating ET in SWAT. The Hargreaves method is simpler than the Penman–Monteith method because it requires fewer input data (maximum and minimum daily temperature) while resulting in a reasonable estimation of ET. The simplicity of the Hargreaves, along with our limited input data, showed the best estimation of ET that led to the SWAT simulated flow closest to historical streamflow conditions. Additionally, as a temperature-based method, the Hargreaves method projects future ET with changes in temperature with CC. After calibration, the simulated SWAT output will be daily streamflow in units of cubic meters per seconds (cms). In addition to examining the statistical distribution of the streamflow metrics we selected, we will also use the Kruskal–Wallis test to see how monthly changes are significantly different between individual and combined modeled scenarios.

To establish a baseline of historical LC scenario, the 2006 (mid-point year within historical streamflow record) National Land Cover Dataset with a 30 m × 30 m cell size was used to run SWAT for this watershed. A total of 18 LC classifications from NLCD 2006 were collapsed into 13 LC classifications in order to be used by the SWAT model for processing (Appendix I). Historic daily climate data were downloaded from the gridMET dataset, a spatially and temporally continuous surface meteorological dataset that are available from 1979 to present (Abatzoglou 2013). Daily maximum and minimum precipitation and cumulative precipitation were the three climate variables extracted and used in the initial run of SWAT. The GridMET dataset has a high spatial resolution of 4 km grids and was validated by Abatzoglou (2013) extensively with surface weather station data. Although the dataset is not good at capturing microclimate under the 4 km spatial scale, it does provide better estimation in SWAT modeling than surface station climate data with more coarse resolution as seen in previous studies (Grusson et al. 2017; Bhattacharya et al. 2019).

After running the SWAT for the historical period of 2002–2016, we divided subsequent SWAT runs into nine different future scenarios for the mid-century 2040–2069 (2050s) and late-century 2070–2099 (2080s) (Table 2). These two periods were selected based on data availability, and the purpose is to break down our analysis further and compare the differences between two periods of climate normal temporally. Future LC datasets were downloaded from the Willamette Water 2100 project website at https://inr.oregonstate.edu/ww2100. Each of the LC datasets in WW2100 was developed specifically based on the three individual GCMs employed in the present paper. These LC scenarios were developed in a continuous feedback modeling system called Willamette Envision that includes both biophysical components and socio-economic components. For example, the vegetation model simulates changes in forest composition, forest area burned by wildfires, and the subsequent impact on timber harvesting and succession of vegetation due to CC (Turner et al. 2015). The WW2100 dataset contains projected LC maps for the period between 2010 and 2100 in the WRB under three different climate scenarios. The three climate scenarios (Table 2) were all part of the CMIP5 GCMs to represent low CC (GFDL-ESM2M), reference CC (MIROC5), and high CC (HadGEM2-ES) (Rupp et al. 2013). With warmer, wetter winter and drier summer, these three future climate scenarios' LCs showed changes in vegetation structure and diversity (Turner et al. 2017). In addition to a gradual increase in developments, upland forest areas are more sensitive under high climate scenarios where mixed forest (+28 to 31%) start to take over a significant portion of evergreen forest (−30 to 44%). In low and reference climate scenarios, mixed forest expands (+4 to 12%) at a slower rate, while evergreen forest declines (−6 to 9%) at a slower rate (Figure 2). A total of six WW2100 LC maps were downloaded and used in SWAT for three climate scenarios in two time periods. LC classification from WW2100 is more explicit in describing different types of vegetation, with expanded categorization in conifer forest age. Reclassification was conducted again to aggregate WW2100 LC class into SWAT land classification for running the model.

We selected the three coupled climate scenarios from a set of 20 CMIP5 models under RCP 8.5 (Table 2). The daily climate data come from the Multivariate Adaptive Constructed Analogs (MACA) dataset, which is a type of the statistical downscaling method for GCMs that uses bias correction procedures and a constructed analogs approach (Abatzoglou & Brown 2012). The downscaled climate data have a spatial resolution of 4 km, capturing near-surface weather conditions for watersheds with complex terrain. The MACA dataset was validated using reanalysis across the western US, showing better performance than linear interpolation-based approaches (Jiang et al. 2018). The selected climate models projected increases in mean annual temperature in the Clackamas river region but has seasonal and inter-model variabilities (Table 3). The models also showed a reduction of seasonal precipitation coupled with a moderate increase in summer temperature, which may induce hydrologic drought, reducing low flows, and threatening summer stream temperature thresholds for aquatic life (Chang et al. 2018). Table 2 illustrates the nine scenarios modeled in our study, with three LCC only scenarios, three CC only scenarios, and three CLC combined scenarios.

Hydrology models often require calibration and validation with observed streamflow data. In our study, we used historical streamflow data between 2002 and 2019 from a USGS gage #14211010 at the mouth of the watershed to calibrate and validate our initial baseline scenario in SWAT. More than a dozen SWAT input parameters were calibrated manually and automatically using SWAT-CUP, a program that is able to run many iterations for a range of parameters and find the best fit values within an acceptable range (Abbaspour 2015). SWAT-CUP has a built-in sequential uncertainty fitting algorithm called SUFI-2. In SUFI-2, uncertainty in input parameters can contribute to output uncertainties; therefore, SWAT-CUP's goal is trying to capture output variables in the 95% probability distributions. SUFI-2 runs several iterations at a time, each time shortening the range of the selected parameters to produce better results in the subsequent iterations. In the SWAT-CUP environment, calibration was done using simulated past daily streamflow against observed daily streamflow data from 2002 to 2014, and the five years, 2014–2019, were then used for validation. Twelve model input parameters chosen based on sensitivity to streamflow and water balance were calibrated, and sensitivity analysis in SUFI-2 was conducted to rank each parameter's weight in the best estimation (Appendix II). We selected three goal functions to evaluate our baseline SWAT model calibration and validation performance. Model efficiency was measured by the coefficient of determination (R²), Nash–Sutcliffe efficiency (NSE) index, and the percent bias (PBIAS).

As shown in Appendix II, model performance was very satisfactory with an NSE value of 0.87 and a R² value of 0.88 after calibration and an NSE value of 0.83 and a R² value of 0.84 after validation on daily streamflow outputs. The PBIAS value was −2.1 and 3.9 for calibration and validation (<25 is considered satisfactory) (Moriasi et al. 2007; Narsimlu et al. 2015), indicating that our model closely estimates the observed daily streamflow with slight underestimation. Our calibration results show slight underestimation during the dry season (June–August) and over-predicts flow during the wet season (November–January), which is a common bias in most hydrologic models (Zhao et al. 2018; Pang et al. 2020). The simulated daily flow allows researchers to capture extreme flow events that are usually within the weekly scale. Parameterization results can be found in Appendix II for the best-fitted value obtained through calibration. The built-in sensitivity analysis algorithm in SWAT-CUP SUFI2 ranks the calibrated parameters based on the response of flow to change in parameter value (Appendix II).

Results showed that under the combined CLC scenario, change in mean monthly flow is higher during the drier season (May–September) and lower in late fall and spring (Figure 3). The upper watershed showed on average slightly higher change (2–6%) in monthly mean flow than the lower watershed. Comparison between two future periods showed that the 2080s have a greater range in the majority of the months. The Kruskal–Wallis test of significance across different individual and combined scenarios showed that the most significant differences are in the wettest months of winter and driest months of summer (Appendix III). However, t-test results showed no statistically significant difference in monthly precipitation between the 2050s and the 2080s (Appendix IV), suggesting uncertainties in choice of GCMs.

LCC only scenarios are expected to decrease by an average of 9% in top 10% flow in the mid- and late-century, while CC and combined scenarios showed an overall pattern of increasing peak flow (2−19%) across all climate scenarios (Figure 4). In combined scenarios, the lower watershed tends to have a slightly greater range (1–2%) in peak flow during the simulated period then the upper watershed. Results of LCC and CC combined scenarios showed similar values to CC only scenarios, suggesting that CC is playing a bigger role in influencing annual peak flow by increasing, while LCC tends to decrease peak flow. Peak flow under low CC scenarios had little changes, while both reference and high CC scenarios showed a wider range in annual top 10% flow values (Figure 4). Between the two simulated future periods, it was surprising to find that the interannual variability 2050s results (s.d. 40.2–74.4) are higher than the 2080s results (s.d. 31.7–64.6). Interannual variability is also generally higher in the lower watershed (s.d. 46.7–74.3) than in the upper watershed (s.d. 30.7–46.7) as well.

Although low flow does not seem to show significant changes across all scenarios, the same pattern of higher variability in the lower watershed (s.d. 5.5–6.2) and lower variability in the upper watershed (s.d. 3.0–6.0) still exist, suggesting the stronger influence of urbanization on low flow. Change in low flow ranges from −20 to +9 cm³/s in the CLC combined scenario. LCC, CC, and combined scenarios all exhibited a decreasing trend in 7-day minimum flow; this is concerning especially during the dry season when water demand is high (Figure 5). In our results, CC only scenarios showed a similar pattern in the magnitude of change as the combined scenarios, once again suggesting that CC is a bigger driver of hydrologic change in the Clackamas river basin. High CLC combined scenarios in the late-century showed the most decrease in 7-day minimum flow. The variability of change across GCMs is highest in the low climate scenario; CC on top of LCC can exacerbate the effect of already decreasing minimum flow.

The mean center-of-volume date for the Clackamas river was found to be shifting earlier across all emission scenarios; both future simulation periods showed a shorter range than historical conditions (Figure 6). Observation of the results showed that in low climate emission scenarios, CT in the 2080s is earlier than the 2050s. However, in reference and high climate emission, this change is reversed with CT shifting later in the 2080s than the 2050s. There are no major differences in the shift between lower and upper parts of the watershed, and CC remains the biggest driver for change in timing of flow in all simulated scenarios.

CC only and combined scenarios showed the most change compared to LCC only scenarios. Streamflow impacts from LCC only scenarios were found to be minimal. Both peak flow and low flow were projected to decrease with LCC, possibly due to the increase of woody wetlands (+12%) and the abundance of mixed forest (+28%). Although the sensitivity of streamflow metrics to LCC was low, peak flow showed somewhat the opposite trajectory of change than CC only scenarios where LCC would decrease peak flow while CC will increase peak flow. These changes in peak flow can be explained by the nature of our study sites, which is heavily dominated by precipitation in the forms of snow and rain, and peak flow events are largely influenced by extreme short-duration precipitation or rain-on-snow events triggered by sudden warming in spring (Safeeq et al. 2015). The CC only scenarios in all streamflow metrics look identical to the combined scenarios, suggesting that CC impacts are driving a big portion of the change even in the CLC combined scenarios. We speculate that in other heavily snow-fed watersheds, basin (large-scale) hydrological change is primarily driven by CC, while local (small-scale) hydrological change is driven by LCC (He et al. 2013; Ahiablame et al. 2017). Highly sensitive streamflow due to climate means the timing, magnitude, and type of precipitation during wet seasons is extremely important for water availability during dry seasons and wet seasons due to concerns for hydrologic extremes (Vano et al. 2015b; Feng & Beighley 2020). Similarly, increasing spring and summer temperatures (2–5 °C) can alter the timing of snowmelt events, causing higher risks of hydrologic drought in dry seasons (Table 3). CC was shown to be impacting runoff on the annual and monthly scale, increasing peak flow, and reducing low flow. The extreme impacts of CC are also dependent on location, where more developed areas can show high sensitivity to streamflow change during extreme events such as winter storms and flooding. Our study confirms that in a highly forested watershed, LCC plays a minimal role compared to CC due to minimal land converted to impervious surfaces (Figure 2). The most significant change observed in this study was the CT of flow. The combined scenarios showing a range of 0–30 days earlier in reaching 50% of annual flow volume can mean that snowmelts will occur much earlier in the year and deplete the already reduced snowpack early enough to cause water availability concerns for managers. The majority of the CRW sits on the colder wetter western cascade that is mountainous and receives rain and snow, while the urbanized part lies on the lower elevation drier Willamette valley. The geologic and geographic differences of the watershed can influence the amount of precipitation and moisture holding properties and can be highly sensitive to climatic changes that will influence the form of precipitation and mechanisms of groundwater recharge. The majority (70%) of our studied watershed lies in western cascade and depends on groundwater during dry seasons to sustain river flow, and it is highly sensitive to climatic changes.

Since we selected the low, mid, and high range of GCMs to represent a different spectrum of scenario outcomes, we expected the high climate scenarios would perform the worst, but that was not the case. Peak flow was more sensitive to reference climate and high climate than low climate. Reference climate becomes more vulnerable to impact in the late-century than in the mid-century. GCM MIROC5 was used for reference climate, and previous research showed that this climate model often underestimates annual precipitation and displayed poor performance at the annual scale with larger uncertainty compared to the HadGEM and GFDL model (Tan et al. 2014). The sensitivity of peak flow to emission scenarios is largely driven by extreme precipitation and snowmelt in our studied watershed caused by an anticipated increase in winter precipitation and earlier snowmelt from warmer spring temperature (Rana & Moradkhani 2016). For low flow index combined scenarios, high climate showed the most decrease in low flow, suggesting a shift in precipitation patterns and highest warming temperature in the summer will have a more drastic effect on low flow during dry seasons than high flow during wet seasons. Variability across climate models is similar, except for a slightly higher range of change in peak flow for the combined scenarios from reference and high climate. As variability increases with high climate scenarios, peak flow becomes more unpredictable during the wet seasons, and flood risk may be greater for the study area and more difficult to forecast and manage. High climate scenarios showed that low flow in the mid- and late-century can be reduced at a magnitude of 0–20 cm s, meaning that low flows will be even lower than current levels in dry seasons, threatening drinking water supply. Drought will be more devastating in the future, and water managers need to seek an alternate source of water, such as creating more reservoirs to meet increasing water demand with urbanization (Moore 2015). Interestingly, there are not many differences in CT change between the climate models, and all three climate models are showing flow volume shifting to earlier in the year. Overall, the reference climate model remains the biggest threat and worst possible outcome in all of the model outputs.

Streamflow is highly sensitive on an LC gradient in watersheds. The lower part of the watersheds with many developments and urban areas showed a very high range of possible outcomes in peak flow and low flow than the upper, forested watershed. This finding matches with previous studies that urbanization can drastically change the local hydrology of the landscape and cause rivers and streams to become more flash during storms (Hale et al. 2016; Anim et al. 2018). Conversion to impervious surfaces can worsen flooding during wet seasons as more precipitation falls as rain than snow (Ansari et al. 2019) in the future as temperature rises. The forested rural subwatershed had the lowest variability across all scenarios modeled, which can be explained by the forecasted LCC where much of the forest vegetation will still exist in the forms of converting from evergreen to mixed forest and woody wetlands. Mixed forest landscapes have the potential to increase runoff and cause streams to be more flashy following precipitation due to the loss of interception and ET from the lack of canopy cover and floor vegetation in fall and winter (Perry & Jones 2017). Woody wetland, on the other hand, can have the benefit of reducing peak flow through increased storage and evaporation. LCC in the upstream section of the watershed will continue to provide essential hydrologic functions like infiltration; however, hydrologic regimes such as annual runoff and peak streamflow may increase from the conversion of conifer to deciduous forest (Tolessa et al. 2017; Qi et al. 2019). The seasonal comparison showed that the sensitivity of flow to the LC gradient is the strongest in late spring and early summer. This result suggests that water managers must allocate water properly to ensure sufficient water available for the peak dry months of July and August.

Running hydrology models have many uncertainties, and scenarios outcomes are highly variable depending on landscape gradients, climate model input, scenario type, and chosen hydrologic models (Hattermann et al. 2018; Chegwidden et al. 2019). Since we conclude that CC is the biggest driver of variability in streamflow indices, we also believe that the three downscaled GCMs selected were the first source of uncertainty in our outputs. We observed mixed results in the forecasted period of the 2050s and the 2080s in streamflow metrics, which is expected based on the common uncertainties that exist in GCMs (Guimberteau et al. 2017). Uncertainty in climate data inputs is likely the cause of reference climate models having a similar performance with high climate models. There are also uncertainties in LC input data as well in SWAT modeling. SWAT is only limited to run one LC map at a time for a given period, and due to the limited temporal resolution of LC projections (1 map per decade), many assumptions were made in running the model assuming static LC for 30 years each time the model was run. Future modeling efforts should incorporate dynamic LC inputs to improve SWAT performance and to reduce uncertainty. Dynamic LC inputs can better capture interannual variability over long periods of simulation; additionally, it can also improve monthly output results (Wang et al. 2018).

Furthermore, compatibility problems with NLCD LC class and SWAT LC class, as well as WW2100 LC classification, required model users to aggregate data and potentially lose essential LC class and original spatial resolution. Lastly, subjective uncertainty could have been introduced during the model parameterization process. Although there are many published studies on instructions to which parameters users should pick to calibrate their studied watershed (Arnold et al. 2012; Narsimlu et al. 2015), many assumptions are made in the calibration process where users choose a uniform value for a certain parameter for the entire watershed. The SWAT model has been widely used on watersheds in many published studies; however, watershed hydrology, geology, and topography can all contribute to model uncertainty. Future improvements to modeling with SWAT should focus on a spatially explicit and dynamic model that can capture the complex interactions between CC and land surface hydrology.

Using the tightly coupled CC and LCC scenarios on a snow-influenced watershed in the PNW over two future periods, we can assess impacts to streamflow indices across a temporal and spatial gradient. Hydrologic variability in both simulated future periods in the CRW is higher in the lower urbanized and agricultural watershed than the upper forested watershed. CC is the primary driver of streamflow change in both future periods across all presumed climate scenarios models. Top 10% flow is projected to increase, while low flow is projected to decrease in combined scenarios across all climate models. CT of flow is expected to shift from mid-May to late-April, due to rising temperature and declining snowpack, threatening water supply in the subsequent dry months. Multiple sources of model limitation and uncertainty include the selection of GCMs, hydrology model parameterization, and LC classification, which are likely causing high variability in streamflow.

Water resources in the CRW are heavily driven by precipitation and snowmelt and are more sensitive to change in climate signals than LC. Our results suggested that future research on similar watersheds should incorporate integrated CC and LCC datasets to better understand the hydrologic sensitivity to each factor. This study showed the need for more adaptive water resource management in the uncertain future and calls for improvements in current modeling techniques to address uncertainties and potential feedback associated with GCMs and LCC.

This research was sponsored by the Institute of Sustainable Solutions at PSU and Clackamas River Water Providers. We appreciate reviewers who clarified many points of the manuscript. Views expressed are our own and do not necessarily reflect those of the sponsoring agencies.

All relevant data are included in the paper or its Supplementary Information.