A geohydrologic framework for characterizing summer streamflow sensitivity to climate warming in the Pacific Northwest, USA

<jats:p>Abstract. Summer streamflows in the Pacific Northwest are largely derived from melting snow and groundwater discharge. As the climate warms, diminishing snowpack and earlier snowmelt will cause reductions in summer streamflow. Most assessments of the impacts of a changing climate to streamflow make use of downscaled temperature and precipitation projections from General Circulation Models (GCMs). Projected climate simulations from these GCMs are often too coarse for planning purposes, as they do not capture smaller scale topographic controls and other important watershed processes. This uncertainty is further amplified when downscaled climate predictions are coupled to macroscale hydrologic models that fail to capture streamflow contributions from deep groundwater. Deep aquifers play an important role in mediating streamflow response to climate change, and groundwater needs to be explicitly incorporated into sensitivity assessments. Here we develop and apply an analytical framework for characterizing summer streamflow sensitivity to a change in the timing and magnitude of recharge in a spatially-explicit fashion. Two patterns emerge from this analysis: first, areas with high streamflow sensitivity also have higher summer streamflows as compared to low sensitivity areas. Second, the level of sensitivity and spatial extent of highly sensitive areas diminishes over time as the summer progresses. Results of this analysis point to a robust, practical, and scalable approach that can help assess risk at the landscape scale, complement the downscaling approach, be applied to any climate scenario of interest, and provide a framework to assist land and water managers adapt to an uncertain and potentially challenging future.</jats:p>

section_content: A fundamental challenge facing scientists and resource managers alike is grounding predictions of climate change and its consequences in specific landscapes and at scales useful for resource planning. This challenge is particularly acute for predictions of water abundance and scarcity, as both the climatic and landscape controls on water availability are typically at a finer scale than representations in the current class of climate and hydrologic models. Resource managers are tasked to plan for an uncertain future by assessing vulnerabilities and sensitivities of different landscapes to change. What strategy should they follow?

One way to assess streamflow vulnerability to changing climate is via a "top-down" approach, which generally involves coupling general circulation models (GCMs) with hydrologic models that predict regional streamflow (e. g., Nash and Gleick, 1991; Hamlet and Lettenmaier, 1999; Nijssen et al., 2001; Christensen et al., 2004; Jha et al., 2004 Jha et al.,, 2006;; Milly et al., 2005; Tohver et al., 2014). This approach has many strengths, including simulation of hydrologic processes under multiple climatic scenarios and across large spatial and temporal scales, and forecasting of hydrographs. 

Published by Copernicus Publications on behalf of the European Geosciences Union. 

But there are also limitations. GCMs coarsely parameterize terrain and fail to incorporate important climatic processes, such as the El Niño-Southern Oscillation and Pacific Decadal Oscillation, in predictions. Higher-resolution regional circulation models (RCMs) that include better topographic representation are improving this situation (Leung and Qian, 2003; Maraun et al., 2010), but accurate forecasts of future climate by this method are still several years off. Moreover, large-scale hydrologic models commonly used in the Pacific Northwest (PNW) for hydrologic forecasting (e. g., variable infiltration capacity -VIC; Liang et al., 1994), do not explicitly simulate streamflow contributions from deep aquifers (Wenger et al., 2010). However, the issue of deep-groundwater representation is not limited to VIC alone. Explicit representation of deep groundwater is approximated by extended soil profiles in many large-scale land surface models (Vano et al., 2012). 

Several recent studies have demonstrated the important role of geologically controlled deep groundwater in mediating streamflow response to climatic variability and warming in the PNW (Jefferson et al., 2008; Tague et al., 2008 Tague et al.,, 2013;; Tague and Grant, 2009; Mayer and Naman, 2011; Waibel et al., 2013). Historical streamflow analysis across the western US underscores the importance of both climatic and geologic controls on streamflow response to climate change (Safeeq et al., 2013). Accordingly, approaches that capture both climatic and geologic controls are needed to identify landscape-level streamflow vulnerability to changing climate. This is particularly critical in the PNW, where local climate, topography and geology combine to dictate hydrologic regimes. 

In the PNW, seasonal asynchrony between winter and spring precipitation and runoff and summer water demand makes summer water supplies scarce and vulnerable (Jaeger et al., 2013). Climate change will intensify this water scarcity by reducing summer streamflows (Safeeq et al., 2013). Declines have the potential to be acute, due to a combination of observed and predicted shifts in precipitation phases from snow to rain, earlier onset and faster rates of snowmelt, and increased summer evapotranspiration (Mote et al., 2005; Stewart et al., 2005; Nolin and Daly, 2006; Das et al., 2011). Increasing interannual variability and changes in extreme flows compound seasonal changes. Luce and Holden (2009), for example, documented widespread declines in the lowest annual flows occurring from 1948 to 2009; these flows are critical for consumptive water use, hydropower, and aquatic biota, including the region's prized and declining salmon populations. 

We present a complementary "bottom-up" approach, focusing on the PNW. Our methodology rests on the analytical framework of Tague and Grant (2009) that characterizes relative summer streamflow sensitivity. Using a rigorous definition of summer streamflow sensitivity as function of the first derivatives of the relationship between discharge and either the timing or magnitude of recharge, we develop a spatial analysis that characterizes summer streamflow sensitivity at a landscape scale. Relationships between observed climate and streamflows at specific gaged locations in diverse hydrogeologic areas are used to extend the sensitivity relationships to ungaged areas and map sensitivity for the entire study region. The uniqueness and strength of this approach is that it is independent of climate change scenarios. Sensitivity is mapped as an intrinsic property of the landscape as interpreted using the average historical climate and other landscape properties, rather than as a response to future climate change alone. 

This sensitivity assessment can then be integrated with climate scenario data to produce regional-scale summer streamflow vulnerability maps. We present an example of how this type of spatial analysis might be applied to National Forest lands in the Pacific Northwest. Land and water managers can tune this type of assessment to their specific needs in order to identify and prioritize actions to adapt to uncertain and potentially challenging future conditions.
section_name: Introduction

section_content: This analysis encompasses Oregon (OR) and Washington (WA) in the northwestern United States (US) with a population of nearly 10. 5 million (US Census Bureau, 2010). The elevation varies from sea level to over 4300 m at Mount Rainier, with the north-south trending mountains of the Cascade Range dividing the western and eastern portions of the states (Fig. 1a ). The study region is divided into 13 physiographic sections (Fig. 1b ) based on common topography, rock type, structure, and geomorphic history (Fenneman and Johnson, 1946). The maritime climate is highly influenced by the Pacific Ocean and varies with elevation and distance from the coast (Fig. 1b and c ). Long-term average precipitation ranges from 150 mm in the Columbia Valley on the east side of the Cascades to ∼ 7000 mm in the Olympic Mountains (Daly et al., 2008, Fig. 1c ). Both OR and WA have extreme wet (winter) and dry (summer) seasons, but the seasonal distribution of precipitation varies between the region's eastern and western halves. While most of the annual precipitation occurs during fall and winter, more frequent summer thunderstorms in the eastern half result in a slightly higher summer precipitation (Mass, 2008 ). An altitudinal temperature gradient, varying by latitude (Fig. 1c ), controls the phase of precipitation with winter rain (R) in lower elevations, seasonal snow at higher elevations (SSZ), and transient snow at intermediate elevations (TSZ) (Jefferson, 2011). The majority of the winter precipitation occurs as rain in the Coast Range and as snow along the Cascades and other ranges (e. g., Wallowa and Blue Mountains). 

This strong climatic gradient and underlying geology that mediate landscape drainage efficiency (Tague and Grant, 2009) are predominant controls on the hydrologic regime of this region (Wigington et al., 2013). For example, streamflow recedes quickly in watersheds with low spring snowmelt and minimal groundwater storage (e. g., the OR Coast Range and western side of the Middle Cascade Mountains known as the Western Cascades) resulting in higher winter peaks and prolonged summer low flows. In contrast, streams in groundwater-dominated regions such as the volcanicdominated central and eastern portion of the Middle Cascade Mountains (High Cascades) show a much more uniform flow regime, with higher summer baseflows, slower recession rates, and significantly lower winter peak flows (Grant, 1997; Tague and Grant, 2004).
section_name: Study location

section_content: Our conceptual model is built around the assumption that the discharge from a watershed depends solely on the amount of aquifer storage. Based on conservation of mass, the water balance within the watershed is given by

where S is water stored in watershed (mm), I R is rainfall (mm day -1 ), I M is snowmelt (mm day -1 ), ET is evapotranspiration (mm day -1 ), Q is discharge (mm day -1 ), and GW in and GW out are the groundwater (mm day -1 ) inflow and outflow, respectively. Change in storage (dS/dt) is positive when

Maximum aquifer storage (dS/dt = 0) occurs when Q = I R + I M + GW in -GW out -ET, which should coincide with peak discharge (dQ/dt = 0) based on the storage-discharge relationship. In reality, since peak discharge always lags the peak recharge (Kirchner, 2009), the peak of the hydrograph will occur when I R + I M + GW in -GW out -ET < Q and thus dS/dt < 0. However, we simplify and assume that at the peak of the hydrograph dS/dt ≈ 0 and hence Q ≈ I R + I M + GW in -GW out -ET, and Eq. ( 1 ) can be simplified to

where Q o is peak discharge (mm). 

The recession curve of the hydrograph, or decay of Q o over time, can be expressed by

where Q(t) is streamflow at time t (in days) from the beginning of the recession period, Q o is streamflow at t = 0, and k is a recession constant (Tallaksen, 1995). As the climate warms, any change in the timing and magnitude of Q o will affect Q(t). Additionally, the recession time t depends on the day of the peak discharge t p and the day t d on which Q is quantified. Hence a more general form of Eq. ( 3 ) can be written as

where Q o and t s are change in peak discharge rate and shift in time driven by climate change, respectively. An earlier shift in peak discharge will result in a negative t s and hence an overall longer recession period between t p and the day t d. Following Tague and Grant (2009), streamflow sensitivities to a shift in magnitude ( Q o ) and timing (t s ) can be described using a first-order derivative of Eq. ( 4 ) with respect to peak discharge Q o and time t:

where terms ε Q o and ε t represent the metrics used in this study to describe the sensitivity of discharge to changes in magnitude of peak discharge and timing, respectively. The negative sign in Eq. ( 6 ) indicates that Q(t) decreases with increasing t. 

The response surfaces of ε Q o and ε t (Fig. 2 ) illustrate the interaction between t and k and how the two sensitivities are expressed over the course of the streamflow recession. In groundwater-dominated systems with low values of k (e. g., High Cascades), ε Q o starts higher at the beginning of the recession and shows a very subtle decline with increasing t (Fig. 2a ). In contrast, in the runoff-dominated systems with high k (e. g., Western Cascades), ε Q o is very comparable to low k systems but diminishes very rapidly with increasing t. In the context of climate change, this suggests that while changes in summer streamflow in groundwaterand runoff-dominated systems with similar t p and Q o may be comparable in the beginning of recession, they vary drastically as the recession progresses. The interaction between t and k for ε t is more complex as compared to ε Q o (Fig. 2b ). In groundwater-dominated systems with low k, ε t starts low and shows a very subtle decline with increasing t. In runoffdominated systems with high k, ε t starts high but diminishes very quickly with increasing t. The very subtle and rapid decline of sensitivities (ε Q o and ε t ) between groundwaterand runoff-dominated systems expressed by the conceptual model are consistent with those expressed in streamflow trends in the empirical record. In groundwater-dominated systems streamflow response to decreasing snowpack is mediated and streamflow continues to decline throughout the summer (Mayer and Naman, 2011; Safeeq et al., 2013). 

Although our conceptual model of streamflow sensitivity is consistent with trends shown in the empirical streamflow record, we recognize that the complexity of the real world is not captured by this simple formulation. Hence, several caveats and assumptions must be emphasized when applying this model. While there is a physical basis for the conceptual model, it is not physically based in a rigorous sense and involves several simplifying assumptions. Watersheds do not typically behave like linear reservoirs; filling (recharge) and emptying (discharge) often occur simultaneously, even during recession periods. Also, groundwater exchange (GW in and GW out ) between watersheds dictates the streamflow regime in some parts of the landscape (Jefferson et al., 2006; Wigington et al., 2013; Patil et al., 2014). This exchange has been excluded from our approach because physically accounting for groundwater gain and loss in this conceptual sensitivity framework with little or no data would undermine its simplicity. This introduces some error in some landscapes, notably those with large groundwater systems in young volcanic terranes. Additionally, this sensitivity approach assumes that Q o and t are independent and any change in Q o will not affect t. This assumption may hold true in rain-dominated systems but could be problematic in snowmelt-driven environments. However, this is much less of an issue in our study domain, where most of the snowmelt occurs during spring and summer recession characteristics depend primarily on peak initial recharge (Q o ). Additionally, approximating the I R or I M for Q o and t R or t M for t p, even when ET ≈ 0 and GW in = GW out (Eq. 2) could result in biased estimates of sensitivity described in Eqs. ( 5 ) and ( 6 ). In places where the reservoir is large, Q o gets delayed following recharge I R or I M, and t R or t M may not represent t p. For example, in watersheds within the seasonal snow zone (see Sect. 4. 2 for the definition), t p is on average delayed by 6 days from t M (Fig. 3 ). In rain-dominated watersheds, the time lag between t R and t p is on average 9 days for the first peak as streamflow recovers from the long summer recession. This time lag between rainfall and streamflow decreases to 1 or 2 days for the subsequent peaks (Fig. 3a, c ). Although the time lag between peak recharge and streamflow may vary significantly depending on data (e. g., observed vs. simulated I M and t M ) and method (e. g., isotopic techniques vs. simple recharge-runoff relationship) used to characterize the relationship, our goal here is to highlight the issue and how it might affect the sensitivity expressed using the conceptual model. Finally, the watershed recession constant, k, may vary year-to-year depending on evapotranspiration losses and other forms of water withdrawals (Thomas et al., 2013), which are not explicitly considered in the model. Given these limitations, our intent is not to precisely predict the change in actual flow regimes, but to assess the comparative sensitivity of those flow regimes across the landscape.
section_name: Conceptual model of streamflow sensitivity

section_name: Parameterizing the model

section_content: Daily average streamflow data for a set of 227 (111 in OR and 116 in WA) unregulated watersheds were obtained from the United States Geologic Survey (USGS) (http:// waterdata. usgs. gov/or/nwis/sw; data accessed on October 31, 2011) and the Oregon Department of Water Resources (http: //apps2. wrd. state. or. us/apps/sw/hydroreport/; data accessed on 1 November 2011) (Fig. 1a ). Watershed drainage areas range from 4 to ∼ 21 000 km 2, with an average of approximately 950 km 2. These watersheds were classified as part of the USGS Hydroclimatic Climatic Data Network (HCDN) (Slack et al., 1993), or were part of the reference gage net-work developed by Falcone et al. (2010) based on Geospatial Attributes of Gages for Evaluating Streamflow (GAGES). Both the HCDN and GAGES data sets have been screened to ensure that they are minimally affected by upstream anthropogenic activities such as irrigation diversions, road networks, and reservoir operations. To minimize the effect of climate bias (i. e., wet vs. dry years) on estimates of k, all selected watersheds were further screened to have a minimum of 20 years of complete daily streamflow data within the water years 1950-2010. Since the majority of the streamflow gages were located in the western half of the study area (Fig. 1a ), we added 12 additional non-reference, non-HCDN gages to the eastern side to ensure a more uniform population of basins. These 12 gages were selected after visual examination of the historic streamflow data records for homogeneity, and review of site information, including a hydrologic disturbance index (Falcone et al., 2010) to ensure there were no major diversions or impoundments. The selected 227 watersheds were delineated using a 30 m resolution digital elevation model (DEM).
section_name: Recession constant (k)

section_content: Following Vogel and Kroll (1992), an automated recession algorithm was employed to search the historical record of daily streamflows for all recession segments lasting 10 days or longer. Peak and end of recession segments were defined as the point when the 3-day moving average streamflow began to recede and rise, respectively. The beginning of the recession (inflection point) was identified following the method of Arnold et al. (1995). To minimize the effect of snowmelt on k, and thereby derive estimates of k that were intrinsic to the geology of the watershed, we excluded recession segments that fell between the onset of snowmelt-derived streamflow pulse and 15 August. We used 15 August as a cutoff for melt-out date determined based on snowpack data from snowpack telemetry (SNOTEL) sites in OR and WA. The date of snowmelt pulse onset was determined following the method of Cayan et al. (2001) and mean flow for calendar days 9-248 after Stewart et al. (2005). Similar to Vogel and Kroll (1992), spurious observations were avoided by only accepting pairs of receding streamflow

where m is the total number of pairs of consecutive daily streamflow, Q t and Q t-1, at each site. Among the 227 watersheds, the values of m varied between 24 and ∼ 8000 (average ∼ 3000). The importance of k in characterizing the low-flow behavior of streams has long been recognized but there is a considerable debate on appropriate techniques for recession analysis (Tallaksen, 1995; Vogel and Kroll, 1996; Smakhtin, 2001; Sujono et al., 2004). Estimates of k are comparable using some techniques (Sujono et al., 2004) but not others (Vogel and Kroll, 1996). To ensure that our k estimates for the candidate sites were robust and not influenced by our choice of the technique for recession analysis, we recalculated k from the master recession curve generated for each site using the matching strip method (Posavec et al., 2006). We also calculated average k from semi-logarithmic plots of individual recession segments lasting 10 days or longer during non-snowmelt periods as described earlier. The recession constant derived from the three methods showed a strong correlation (R > 0. 77, p < 0. 001). We used the recession constant k from Eq. ( 7 ) in the sensitivity analysis.
section_name: Recession analysis

section_content: We established a regression model for transferring k to the ungaged landscape. Average watershed relief and slope were estimated from a 30 m DEM using the ArcGIS spatial analyst. Soil permeability (K soil, cm h -1 ) values for the top 10 cm soil depth were obtained from the STATSGO database (Miller and White, 1998 ; available online: http://www. cei. psu. edu). A digital 1 : 500 000 scale ArcGIS coverage of aquifer permeability (K aqu, m day -1 ) derived from existing aquifer unit maps for eastern OR (Gonthier, 1984) and western OR (McFarland, 1983) was obtained from Wigington et al. (2013). Because this K aqu data set was not available for WA, we developed a geologic index (ranging from 1 to 9, with higher values corresponding to higher permeability) for OR and WA based on a 1 : 500 000 scale aquifer porosity and rock unit map (Huntting et al., 1961; Walker et al., 2003). A regression between drainage densities estimated using the National Hydrography Dataset (NHD) flowlines and the area-weighted geologic index was used to assign the K aqu values to each geologic index in WA. Area-weighted values of average relief, slope, K soil, and K aqu were determined and log-transformed prior to the regression analysis. 

Starting with the entire list of parameters (i. e., relief, slope, K soil, and K aqu ) from 227 watersheds, we developed a multiple linear regression model. The established regression model was then used to generalize k values across the region (wall-to-wall) at the landscape scale. The prediction for k was made at the fifth field Hydrologic Unit Code (HUC) scale of the national Watershed Boundary Dataset; fifth field HUC units are termed watersheds and typically range in area from 160 to 1010 km 2. Outliers in the model parameters were identified based on Cook's distance (Cook, 2000) and subsequently excluded from the regression analysis using the recommended threshold of 4/n sn i -1, where n s is the sample size and n i is the number of independent variables. Non-significant (p ≥ 0. 15) model parameters were then eliminated via backward stepwise regression, until all remaining parameters were significant and the predictive power of the equation (based on adjusted R 2 ) began to decline. This regression equation was developed individually for OR and WA as well as the entire domain with both states combined (Table 1 ). The correlation matrix for the watershed parameters used for predicting k showed strong cross-correlation (as high as 0. 72), particularly among K aqu, Slope, and K soil in OR. However, since these variables are used to predict k and not to characterize their relationship with each other, the cross-correlation and sign of the regression coefficients can be ignored. 

The regression coefficients (R 2 ) for the three geographic domains (OR: model 1a, WA: model 1b, or OR and WA combined: model 2) ranges between 0. 44 for WA and 0. 59 for OR (Table 1 ), which are within the range of values reported elsewhere with a different set of independent variables (e. g., Thomas et al., 2013). The overall standard error of the estimate is low for the fitted regressions, and modeled k is only slightly biased, overpredicting small values and underpredicting higher values of k (Fig. 4 ). There is no clear spatial pattern of systematic bias based on residuals, however (Fig. 4d ). The predicted k map using model 2 at the fifth field HUC scale broadly distinguishes among different hydrologic regions with different drainage characteristics, including fast-draining regions such as the OR Coast Range, parts of the Columbia River basin in OR and WA and the Owyhee uplands and much of the Ochoco Mountains in OR. Slower-draining regions include the eastern (High) Cascades in OR and WA and the Okanogan highlands in WA (Fig. 5a ), but the Okanogan k values are at the high end of the range for this bin (0. 02-0. 04).
section_name: Regression model development

section_content: We approximated the peak discharge (Q o ) in Eq. ( 2 ) by peak recharge (I R or I M depending on the dominant recharge type), ignoring the groundwater exchange between HUC units (GW in = GW out ≈ 0) and with ET ≈ 0 at the start of the recession. In the PNW, the peak recharge pulse during the water year can be either rain or snowmelt, depending on geographic location. We assigned the primary type of peak recharge pulse (rain or snowmelt) based on a temperature threshold and snow to precipitation proportion. Following Jefferson (2011) and Nolin and Daly (2006), a winter temperature-based threshold of 0 • C was chosen to approximate the boundary between the transitional snow zone (TSZ) and rain zone, while -2 • C was chosen to approximate the boundary between the seasonal snow zone (SSZ) and TSZ. Following Knowles et al. (2006), we define winter as beginning in November, rather than January, and only use wet-day minimum temperatures, which showed a strong correlation with the snow to precipitation ratio. We defined a wet-day as a day when daily precipitation is greater than zero. In addition, we used the temperature threshold-based empirical relationship of Dai (2008) and the United States Army Corps of Engineers (USACE, 1956) to calculate the median value (water year 1916-2006) of the fraction of annual precipitation falling as snow. We classified the peak recharge pulse as rain for the entire area within the identified rain zone and the portion of area in TSZ with a median snow fraction < 10 %; the remaining TSZ and entire SSZ were classified as snowmelt recharge pulse (Fig. 5b ). A lack of spatially distributed precipitation gauge and snowpack telemetry sites, particularly at higher altitudes, precluded using empirical data to calculate recharge magnitude and timing. Instead, we calculated the peak recharge magnitude (I R and I M ) and timing (t R and t M ) using spatially distributed gridded (1/16th degree resolution) daily precipitation and VIC-simulated daily snowmelt data from Hamlet et al. (2013). The simulated snowmelt data from Hamlet et al. (2013) were limited to the Columbia River basin and coastal river basins of OR and WA and did not include the OR portions of the Klamath and Great basins. VIC-simulated daily snowmelt data for the Klamath and Great Basins at 1/8 • spatial resolution were obtained from the US Bureau of Reclamation (Reclamation, 2011). VIC uses a two-layer energy and mass balance approach to model the process of snow accumulation and melt; descriptions of snow accumulation and melt processes within the VIC model are well described elsewhere (Liang et al., 1994; Ni-Meister and Gao, 2011). 

The daily (1-365) average ) maximum 1-day recharge, I R and I M were calculated on the water year basis as

where R is the daily precipitation (mm), M is the daily snowmelt (mm), and N is the length of record (year). The corresponding timing t R and t M were calculated as the day of water year on which I R and I M occurred. The spatial distribution of recharge magnitude (I R and I M ) and timing (t R and t M ) show distinct geographic contrasts between the eastern and western study domains (Fig. 6 ). The average peak daily recharge from precipitation (I R ) varies from less than 5 mm day -1 in the Walla Walla Plateau and much of eastern OR to as high as 44 mm day -1 in the Olympic Mountains to the west. Similarly, the average daily peak snowmelt (I M ) varies between 0 in coastal southeastern OR to as much 40 mm day -1 in northern WA. Although the magnitudes of I R and I M are small in northeastern WA and much of eastern OR as compared to those in the Coast Range, northern WA, and Cascades, they occur later during the water year. In northern WA, the timing of I M occurs quite late during the water year (Fig. 6 ). Timing of I R is also quite variable across the region and occurs as early as October to as late as mid-September (Fig. 6 ). For the sensitivity analysis, in systems with rain as dominant recharge we substituted Q o with I R and t p with t R. Similarly, in systems with snowmelt as dominant recharge we substituted Q o with I M and t p with t M.
section_name: Historical recharge magnitude and timing (Q o , t p )

section_content: Changes in actual streamflow in the future will not only depend on the intrinsic sensitivity of the landscape but also the magnitude and direction of climate change in terms of magnitude (I R or I M ) and timing (t R or t M ) of recharge to which a landscape is exposed. The actual exposure or magnitude of change in I R or I M and t R or t M will depend on future emission scenarios, which are highly uncertain. However, to illustrate this concept of intrinsic sensitivity and exposure, we present a climate change scenario consistent with regionalscale climate projections for the PNW of decreasing snowpacks (Mote, 2003; Elsner et al., 2010) as a proxy for exposure. An integrated daily snow product based on the 1 km resolution Snow Data Assimilation System (Carroll et al., 2001) was selected and I M and t M were calculated as described earlier. We used the differences between I M and t M values for the wet year 2004 (an El Niño year) and the dry year 2011 (a La Niña year), which correspond to a ∼ 50 % regional snowpack decline, as a potential climate change scenario. Changes in precipitation magnitude and timing are unclear for this region (Salathe et al., 2007; Mote and Salathe, 2010), and were excluded from this analysis.
section_name: Future recharge magnitude and timing (Q o , t p )

section_content: We validated our derived streamflow sensitivities (ε Q o and ε t ) against empirical measures of climate sensitivity extracted from historical records of 217 (Fig. 1a ) watersheds for the months of July, August, and September. Our approach was to use streamflow response to historical climate extremes as a proxy for streamflow sensitivity. Measures used included the (1) change in streamflow with respect to a change in annual precipitation between wet and dry periods; and (2) change in streamflow with respect to a change in spring air temperature between cool and warm periods. These two empirical measures of sensitivity were calculated as

Average annual precipitation (P ) for each watershed was used to identify the 5 years with the lowest and highest precipitation as dry and wet periods, respectively. Similarly, the watershed average of mean daily spring (April-June) temperature (T ) was used to identify the 5 years with the coolest and warmest springs. This approach is analogous to the precipitation and temperature elasticity measure of streamflow sensitivity proposed by Schaake (1990) and Sankarasubramanian et al. (2001). The empirical measures ε p and ε T were calculated as an indicator of streamflow sensitivity to a change in magnitude and timing of recharge, respectively. However, magnitude (I R and I M ) and timing (t R and t M ) are each affected by wet and dry periods and cool and warm springs (Table 2 ). Also, the effect of wet and dry climate on peak recharge magnitude and timing differs for rain-and snowmelt-dominated systems. For example, during a wet as compared to a dry period, t M shifts 16 days later, whereas t R shifts 20 days earlier. Hence, the empirical measures ε p and ε T are representative of the streamflow sensitivities as a convolution of timing and magnitude. We used the nonparametric Spearman rank correlation (ρ) coefficient to evaluate the correspondence between empirical (ε p and ε T ) and conceptual (ε Q o and ε t ) measures of streamflow sensitivities. Spearman rank correlation is less sensitive to outliers and considered a robust alternative to the Pearson product moment correlation. 

6 Results and discussion
section_name: Model validation

section_content: Summer streamflow sensitivities derived from the conceptual framework are in agreement with the climate sensitivity estimators calculated from historical data (Table 2 ). The absolute magnitudes of both empirical (ε p and ε T ) and conceptual (ε Q o and ε t ) measures of streamflow sensitivities decrease from July to September. Also, both precipitationand temperature-based estimators of streamflow sensitivity ε p and ε T are significantly (p < 0. 001) correlated with ε Q o and ε t. The Spearman rank correlation coefficient for ε p and ε Q o decreases from 0. 73 in July to 0. 50 in September, and for ε p and ε t decreases from 0. 77 in July to 0. 54 in September. 

The Spearman rank correlations between ε p and ε Q o or ε t are weaker and ranged between -0. 66 (ε p vs. ε Q o ) and -0. 71 (ε p vs. ε t ) in July and -0. 5 in September. The overall slightly lower values of Spearman rank correlations between empirical and conceptual measures of streamflow sensitivities are not surprising given the fact that changes in I M and t M between wet and dry periods were very small. Similarly, between cool and warm periods I R and t R were relatively constant. So although we used a total of 217 watersheds for validation, not all of them were subjected to a change in magnitude and timing of recharge between wet and dry or cool and warm periods. In fact, all of the rain-dominated watersheds had similar I R and t R between cool and warm periods. This smaller change in I R and t R limits the range of our validation for rain-dominated watersheds.
section_name: Sensitivity validation

section_content: Streamflow sensitivities to a change in magnitude, ε Q o, are very similar during the first weeks after peak recharge for all HUC units across the range of k values (Fig. 7a ). In groundwater-dominated HUCs, the ε Q o are mediated and show very sharp contrasts from runoff-dominated HUCs, even after 110 days of recession. Since peak recharge I M occurs late during the year in most of the low k HUCs (Fig. 6 ), these mediated sensitivities will be expressed throughout the summer. In contrast, the sensitivities to a change in timing, ε t, are very different during the first weeks after peak recharge across all HUC units (Fig. 7b ). Most of the HUCs with higher ε t (> 0. 5 mm day -1 ) are in the rain-dominated Coast Range (Fig. 1 ), where recharge magnitude (I R ) is higher overall when compared to the snow-dominated Cascades, Olympics, and other western parts of OR and WA. However, in most of these coastal HUCs the peak recharge occurs early in the year (Fig. 6 ), resulting in a long recession with lower sensitivities in the summer months. Summer streamflow sensitivities to a change in the magnitude (ε Q o ) and timing (ε t ) of recharge at the beginning of July, August, and September show several distinct patterns (Fig. 8 ). First, there is a clear north-south grain to the sensitivity of both variables due primarily to the corresponding orientation of the topography, with the Cascade Range in both OR and WA clearly showing up as most sensitive to both types of changes. Snow-dominated regions with late melt, such as the mountains along the WA-Canada border and the Wallowa Mountains in OR also show a high, though diminished, sensitivity. Second, the maps show that fifth field HUCs sensitive to a change in magnitude (I R and I M ) are also sensitive to a change in timing (t R and t M ). Third, the level of sensitivity and its spatial extent diminish as the day of interest (t d ) moves from early to late summer. The highest magnitudes of sensitivity to changes in I R and I M, were 0. 47, 0. 25, and 0. 14 mm mm -1 at the start of July, August, and September, respectively; The highest magnitudes of sensitivity to changes in t R and t M were 0. 28, 0. 10, and 0. 03 mm day -1, at the start of July, August, and September, respectively. The highest sensitivity for July streamflow is primarily located in the northern WA and along the Cascades, but portions of OR Cascades continue to show high sensitivity throughout the summer. This contrasting pattern is attributed to relatively high k values in the OR Cascades compared to northern WA. By the end of August, OR Cascade streams are mainly sourced from deep groundwater, as most of the aboveground storage in the form of snow has melted out (Tague and Grant, 2004). 

The influence of k becomes more important than peak recharge magnitude and timing as summer proceeds. Thus, although the different regions display similar levels of sensitivity, the reasons for this sensitivity vary by locale. In contrast, summer streamflow (i. e., July, August, and September) in HUCs that receive recharge in the form of rain (e. g., Coast Range) and do not have deep groundwater are less sensitive to a change in the I R or t R compared to HUCs driven by snowmelt recharge (e. g., High Cascade range and much of northern WA). This lower sensitivity primarily results from peak rainfall occurring earlier in the year (Fig. 6 ), leading to a long summer recession. A similar low sensitivity is observed in eastern OR, where peak snowmelt occurs later in the year, but the magnitude of recharge I M is small and there is very little deep-groundwater contribution to sustain the recession. Over the entire study area, streamflow at the start of July is at least moderately sensitive (ε Q o and ε t > 0. 001) to a change in peak recharge magnitude and timing in 49 and 27 % of the area, respectively. As the day of interest moves towards the start of September, the spatial extent of at least moderately sensitive areas diminishes to 25 and 11 % of the region for ε Q o and ε t, respectively. Within the individual states, streamflow at the start of July in OR is at least moderately sensitive in 38 and 16 % of the area, as compared to 64 and 44 % of the area in WA, to a change in peak recharge magnitude and timing, respectively. Similarly, streamflow at the start of September in OR is at least moderately sensitive in 15 and 6 % of the area as compared to 39 and 18 % of the area in WA, to a change in peak recharge magnitude and timing, respectively.
section_name: Sensitivity analysis and distribution

section_content: This analysis yields a spatially explicit prediction of the sensitivity of late summer streamflow to climate change based on the convolution of geology, as represented by k, and recharge dynamics, as represented by I R, I M, t R and t M (Fig. 8 ). To better understand this sensitivity, we consider how the processes driving it vary across the landscape. For example, the OR High Cascades and much of WA show similar levels of sensitivity, but for different reasons. The OR High Cascades are sensitive because of low k and, as a result, abundant deep, and slow-draining groundwater that recharges streams over many months. Peak snowmelt recharge, I M in much of the OR Cascades, is not only small compared to northern WA, but also melts earlier (Fig. 6 ), leaving deep groundwater as the only source of late season streamflow. These groundwater-dominated landscapes in effect "remember" changes in climate as reflected in either the magnitude or timing of recharge in the winter or spring, resulting in higher sensitivity of late-season streamflow. 

In contrast, much of northern WA is sensitive not because of low k but because of higher I R or I M and late t R and t M. The I M is higher in much of this region and melts later during the year (Fig. 6 ), contributing a substantial portion of the late season streamflow. If the climate changes so that less snow accumulates and snowmelt occurs earlier in spring, the corresponding changes in recharge timing and magnitude are reflected in late summer streamflow, which relies almost exclusively on snowmelt in this region. 

The hydrogeologic sensitivities (Fig. 8 ) illustrate the magnitude of change to existing summer streamflows during early July, August, and September, per unit change in recharge magnitude and timing. Hence, the sensitivity is an intrinsic, mappable landscape property driven primarily by current climate and geology. This information is valuable for climate change planning and mitigation efforts, particularly in ungauged basins, which represent most of the landscape. Our analysis predicts sensitivity to change, but not actual changes to magnitude or timing of streamflow. Actual changes in summer streamflow are a product of both this hydrogeologic sensitivity (Fig. 8 ) and realized changes in I R or I M and t R or t M under a given climate change scenario. Changes in ET are also a factor, but are not considered here. 

Summer streamflow change resulting from this test scenario can be expressed both in absolute (units of flow increase or decrease over time) and relative (percentage increase or decrease over time with respect to Q o ) terms, depending on the application and subject of interest. The average change in I M and t M between the year 2004 and 2011 was 4. 1 ± 4. 5 mm and 38 ± 34 days, respectively. We then calculated late summer streamflow at the beginning of July, August, and September using the change in I M and t M values separately (Fig. 9 ). Only 7 % of the region showed a decline in 1 July streamflow by at least 1 mm (a threshold equivalent to average daily September streamflow) under the I M scenario as compared to 8 % under the t M scenario. Most of the HUCs with a 1 mm or greater decline are located in WA. Nearly 16 % of the area in WA showed at least a 1 mm decline in 1 July streamflow, as compared to only 3 % in OR, to a change in t M between the years 2004 and 2011. Similarly, 12 % of the area in WA showed at least a 1 mm decline in 1 July streamflow as compared to only 3 % in OR to a change in I M between the years 2004 and 2011. As expected, streamflow changes in July were larger than in August and September under both the I M (Fig. 9a ) and t M (Fig. 9b ) scenarios. Relative changes (%) in streamflow were calculated after normalizing the absolute change by the peak snowmelt recharge (I M ) as a proxy for Q o. In the absence of spatially distributed observed streamflow data, we utilized the peak recharge as a proxy for available water in the streams at the start of the recession. In general, areas showing greater absolute change also showed greater relative change (Fig. 9a and b ). 

This disparity between absolute and relative change across the landscape illustrates a key aspect of interpreting sensitivity: our prediction of future streamflows reflects both the intrinsic sensitivity of the landscape (as reflected in k and average historic climate) as well as changes in snowpack between cooler and warmer years. Both factors affect the timing or magnitude of recharge. Specifically, under our assumed scenario, the changes in I M and t M are greater in places with "warmer" snowpacks (Nolin and Daly, 2006), such as the Cascades and other mountain ranges that are closer to marine influence (e. g., Olympics, Fig. 1b ). In these areas, small temperature changes directly affect the total proportion of snow to precipitation. In contrast, colder snowpack areas such as the Harney and Great Basins, Payette, and Walla Walla Plateau (Fig. 1b ) are less sensitive to tem-perature changes. The net effect to streamflow is that some regions (e. g., Northern Cascades, Fig. 1b ) experience both more vulnerable snowpack and more sensitive landscapes (i. e., lower k values). This is reflected in both a greater absolute and relative change (Fig. 9 ). The drier eastern portions of the study region, in contrast, have lower absolute change because their snowpacks are relatively insensitive to warming, and k values are higher.
section_name: Summer streamflow vulnerability

section_content: A central goal in developing this spatially explicit, analytical framework was to help resource managers, such as the US Forest Service (USFS), evaluate vulnerabilities of key resources to changing summer streamflows, and develop and implement adaption strategies to reduce potential impacts. While such strategies may introduce some new activities (e. g., facilitated migration of species, mulching forests) (Grant et al., 2013), we expect that most will involve adjustments in the location, timing, and scope of current actions or modification of their site-specific designs. 

To explore this, we consider how this type of spatial analysis might inform management of National Forest lands in the Pacific Northwest. National Forests comprise a particularly 1 ) the intrinsic sensitivities to changes in peak snowmelt magnitude (Fig. 8 ); and (2) a scenario similar to the differences experienced between a warm, dry year (2003, El Niño) and a cool, wet year (2011, La Niña). Gray areas are rain-dominated recharge and were excluded from this analysis. (b) Predicted decline in streamflow in absolute (i) and relative (ii) terms, based on (1) the intrinsic sensitivities to changes in peak snowmelt timing (Fig. 8 ); and (2) a scenario similar to the difference experienced between a warm, dry year (2003, El Niño) and a cool, wet year (2011, La Niña). Gray areas are rain-dominated recharge and were excluded from this analysis. large fraction of the region (nearly 27 % of OR and WA) and support diverse, valuable, and climate-sensitive resources. The largest changes in summer streamflows are expected to occur on these forest lands, which may affect and alter numerous forest management activities. Such activities include timber harvest and fuels management, watershed restoration, assessment and monitoring, and construction and operation of dams, water diversions, roads, and recreational facilities. 

Watershed restoration is currently a major focus for the USFS (Potyondy and Geier, 2011). Much of this work in the Pacific Northwest is directed towards maintaining or improving water quality and aquatic habitats for salmon and other cold water biota, as directed by the Northwest Forest Plan and other forest plans in the region. Common restoration actions include removal of physical barriers in streams (e. g., poorly designed culverts), road improvements and decommissioning, improved livestock management, reconstruction of stream channels and floodplains, restoration of riparian vegetation and streamflows, decommissioning or alteration of dams and water diversions, and enhancement of instream habitats via additions of wood, boulders, and nutrients (Roni et al., 2002). 

Implementing these restoration projects in a "climateinformed" way is critical, as changes in summer streamflows and other habitat components (e. g., stream thermal regimes) may significantly influence their effectiveness (Battin et al., 2007). This can be accomplished by integrating assessment products like the one presented here into existing strategic planning and project design processes. For example, to maximize the effectiveness of its restoration program, the USFS is currently focusing investments in "priority watersheds" based on assessments of non-climatic stressors and other factors (Watershed Condition Framework at http://www. fs. fed. us/publications/watershed/). In the PNW, those watersheds where the greatest ecological gains can be achieved with the least funding have typically been selected as priorities. In general, such areas have high ecological values (e. g., high biodiversity, rare or legally protected species), mild to modest levels of non-climatic impacts (e. g., water diversions, water quality problems, altered stream habitats), high sensitivities to those impacts (e. g., cold water biota with narrow thermal tolerances), and significant opportunities for restoration (e. g., important and technically solvable problems, sufficient financial resources and workforce capacity, community support, few legal barriers). 

This sensitivity assessment provides an opportunity to consider an additional factor in the priority-setting: climateinduced changes in summer streamflow. In many cases, such changes may not alter priority areas selected for restoration. For example, current priority watersheds may remain priorities after consideration of climate change information (Fig. 10 ). In others, however, likely climate impacts may shift emphasis away from some watersheds and towards others. For example, watersheds with large projected changes in summer streamflows and water resources highly sensitive to those changes may be considered a lower restoration priority if restoration treatments are unlikely to address the cumulative effects of both climatic and non-climatic impacts or if the cost of those treatments greatly exceed available funding (i. e., adaptive capacity is limited). Conversely, the relative priority of other watersheds may increase in cases where significant climate impacts are expected, but managing both climatic and non-climatic impacts is deemed technically, socially, and financially achievable (Fig. 10 ). 

Moreover, this analysis could influence the type, intensity, location, or timing of restoration actions considered necessary to sustain critical resources in priority watersheds, both at a watershed-and project scale. The prospect of late-season streamflow change in some portions of the watershed could lead to redesign of water diversions, proactive efforts to reduce stream temperatures, rethinking low-flow channel dimensions for fish passage and stream channel reconstruction projects, and reconsideration of what riparian species are likely to survive into the future.
section_name: Management applications

section_content: Our results provide a hydrogeologic framework to identify watersheds most and least vulnerable to summer streamflow changes. This method reveals landscape -level patterns and their relationship to topographic, geologic and climatic controls, and can be applied to interpret the effects of any climate change scenario of interest. As such, we believe the sensitivity maps represent a robust, scalable tool that can be used in climate change assessment and adaptation in both gaged and ungauged basins. 

Lack of geologic (i. e., aquifer permeability) and snowmelt information at appropriate spatial scales and accuracies to predict drainage efficiency and peak recharge magnitude and timing is a challenge. For example, aquifer permeability used for OR and WA at the scale of 1 : 500 000 reflects limited spatial heterogeneity and it is unclear how a finer-scale (i. e., 1 : 100 000) permeability or geology map would influence k. Similarly, we relied on simulated historic snowmelt data at 1/16 and 1/8 • grid resolution due to the absence of long-term, spatially distributed measurements. It is unclear how the changes in temperature and precipitation will affect our assumption to approximate the peak discharge with recharge. As the climate continues to warm, the time lag between recharge and streamflow (Fig. 3 ) in rain-and snow-dominated watersheds will likely shift. More rain instead of snow will also alter the dominant recharge regime (Fig. 5b ) and eventually the streamflow sensitivities. Also, this sensitivity analysis should be applied carefully in places where subsurface groundwater exchange or summer evapotranspiration dominate summer streamflow regime. As finerresolution data on both geological and climatic factors becomes available, this approach can be refined to capture new information. More broadly, we recognize that this approach does not yield the specific streamflow values or future hydrographs of the current generation of hydrologic models. There are many applications where having a spatio-temporal prediction of how much water is present would be quite useful. Beyond the uncertainty in both our approach and streamflow modeling, each method has strengths and limitations. The spatial map of sensitivity reveals broad landscape patterns and is applicable where data, time, or cost limit applying a more sophisticated hydrologic model. Hydrologic models give detailed predictions, but may not always illuminate underlying mechanisms or provide sound future predictions. Both approaches have their place. Although our results are independent of GCM predictions, the two approaches are not necessarily mutually exclusive. New CMIP5 high-resolution, terrain-sensitive model predictions could be incorporated into this framework. 

Predicting future streamflows is an uncertain task at best, but is essential to address a rapidly changing environment. The "bottom-up" approach described here is intended to complement other "top-down" approaches involving sophisticated and coupled climate and hydrologic models. These spatial maps, based on simple theory and supported by empirical data, represent spatially explicit hypotheses about how streamflow is expected to respond to climate changes in the future. Other, more complex approaches also yield spatially explicit hypotheses in the form of future hydrographs. We can now compare these two approaches, highlight their strengths and limitations, and integrate knowledge from each to guide managers and communities in facing the uncertain future of water resources in the Pacific Northwest and beyond.
section_name: Conclusions

section_content: Acknowledgements. The authors gratefully thank three anonymous reviewers for their helpful comments. They acknowledge funding support from the Oregon Watershed Enhancement Board, the Bureau of Land Management (Oregon) and the USDA Forest Service Region 6 and Pacific Northwest Research Station. 

Edited by: M. Mikos