Climate change impacts on snow and streamflow drought regimes in four ecoregions of British Columbia

Abstract In many regions with seasonal snow cover, summer streamflow is primarily sustained by groundwater that is recharged during the snowmelt period. Therefore, below-normal snowpack (snow drought) may lead to below-normal summer streamflow (streamflow drought). Summer streamflow is important for supplying human needs and sustaining ecosystems. Climate change impacts on snow have been widely studied, but the relationship between snow drought and streamflow drought is not well understood. In this study, a combined investigation of climate change impacts on snow drought and streamflow drought was completed using generic groundwater – surface water models for four headwater catchments in different ecoregions of British Columbia. Results show that, in response to increased precipitation and temperature, the snow drought regime changes substantially for all four catchments. Warm snow droughts, which are caused by above-normal winter temperatures, increase in frequency, and dry snow droughts, which are caused by below-normal winter precipitation, decrease in frequency. The shift toward more frequent and severe temperature-related snow droughts leads to decreased summer streamflow, decreased summer groundwater storage, and longer, more severe summer low flow periods. Moreover, snow droughts propagate into summer streamflow droughts more frequently in the future time periods (2050s, 2080s) as compared to the baseline 1980s period. Thus, warm snow droughts not only become more frequent and severe in the future but also more likely to result in summer streamflow drought conditions.


Introduction
If temperatures rise as expected (Intergovernmental Panel on Climate Change (IPCC), 2013), precipitation will be more likely to fall as rain than as snow. Warming alone is expected to have large impacts on the hydrologic regimes of catchments with seasonal snow cover (Barnett, Adam, and Lettenmaier 2005), with decreased annual snowpack leading to earlier snowmelt and diminished and potentially warmer late summer flows (Barnett et al. 2008;Seager et al. 2013;Godsey, Kirchner, and Tague 2014;Reynolds, Shafroth, and Poff 2015;Service 2015;Jenicek et al. 2016). Snowmelt is more effective than rainwater at infiltrating beyond the root zone (Earman et al. 2006), and the spring snowmelt pulse often makes up a large fraction of the groundwater recharge in seasonally snow-covered catchments (Winograd, Riggs, and Coplen 1998;Earman et al. 2006;Ajami et al. 2012). This pulse of groundwater recharge from snowmelt is an important component of the hydrologic regime, and the spring groundwater recharge plays a key role in sustaining streamflow during the summer Tague et al. 2008;Tague and Grant 2009;Safeeq et al. 2013;Kormos et al. 2016;Cooper et al. 2018).
A lack of snow accumulation in winter, i.e. a 'snow drought', thus has important implications for water quantity in the following summer season when water demands are high. Snow drought (Ludlum 1978;Wiesnet 1981) can be caused by below-average precipitation (dry snow drought), above-average winter temperatures (warm snow drought), or both (warm and dry snow drought) (Harpold, Dettinger, and Rajagopal 2017). All snow drought types can lead to below-average summer streamflow (Harpold, Dettinger, and Rajagopal 2017;Dierauer, Whitfield, and Allen 2018). Streamflow droughts can propagate directly from snow droughts, with warm and dry winters producing longer, more severe low flow periods in the summer (Dierauer, Whitfield, and Allen 2018). The role of temperature in snowmelt hydrology has been widely studied (Leith and Whitfield 1998;Whitfield and Cannon 2000;Adam, Hamlet, and Lettenmaier 2009;D ery et al. 2009;Pederson et al. 2011;among others), and many studies have documented the importance of snow for groundwater recharge and/or summer streamflow (Green et al. 2011;Huntington and Niswonger 2012;Godsey, Kirchner, and Tague 2014;Markovich, Maxwell, and Fogg 2016;Evans et al. 2018). No studies have explicitly quantified the impact of different snow drought types, i.e. dry, warm, or warm and dry, on the severity and duration of summer streamflow droughts.
Recent work by Dierauer, Allen, and Whitfield (2019) showed that temperature-related snow drought risk increases with increasing mean winter temperatures and identified temperature thresholds above which hydroclimatic change 'accelerates'. As temperatures rise, increased frequency and severity of temperature-related snow drought will likely lead to increased frequency and severity of summer streamflow droughts (Harpold, Dettinger, and Rajagopal 2017;Dierauer, Whitfield, and Allen 2018). This study aims to test the hypothesis that the magnitude of these changes will depend on a catchment's starting point, i.e. its baseline mean winter temperature (Nolin and Daly 2006;Tennant et al. 2015). The specific objectives of this study are as follows: 1. Use generic groundwater-surface water (GW-SW) models to quantify the impact of climate change on (a) the annual and intra-annual water balance, (b) the snow drought regime, and (c) the seasonal (summer versus winter) low flow magnitude for catchments that span a large range of baseline climate conditions. 2. Determine if different snow drought types (dry, warm, or warm and dry) have significant impacts on summer and winter low flows. 3. Quantify the historical and projected frequency of hydrological drought types that propagate directly from snow droughts and determine if changes in the drought regime relate to baseline climate conditions.

Study locations
Four headwater catchments spanning a range of climate conditions were chosen for this study. Each is located in a different level I ecoregion, a broad region of general similarity in environmental resources (Commission for Environmental Cooperation (CEC) 2011) of British Columbia (Table 1 and Figure 1), and each represents a municipal, agricultural, or industrial surface water supply source. The Fort Nelson River headwater catchment is the coldest and has a mean annual temperature of À0.6 C. The Fort Nelson catchment is located in the Hay and Slave River Lowlands of the Taiga ecoregion. The Blueberry River headwater catchment is the second coldest and has a mean annual temperature of 0 C. The Blueberry catchment is located in the Clear Hills and Western Alberta Uplands of the Northern Forests ecoregion. Both catchments have relatively cold, dry winters, and, on average, receive less than 15 cm of snow per year. Additionally, both catchments are in Northeast British Columbia (NEBC)an area of intense shale gas development where multi-stage hydraulic fracturing operations require large quantities of water. Whiteman Creek headwater catchment has a mean annual temperature of 2.1 C and is located in the Thompson-Okanagan Plateau of the North American Deserts ecoregion. The Thompson-Okanagan Plateau has a dry continental climate and is in the rain-shadow of the Coast and Cascade Mountain Ranges; however, the region supports a strong agricultural industry that has a high irrigation demand (Neilsen et al. 2006). The Capilano River headwater catchment is the warmest catchment in this study and has a mean annual temperature of 5.9 C. It is located in the Pacific and Nass Ranges of the Marine West Coast Forests ecoregion. The Pacific and Nass Ranges have a wet maritime climate, and the headwater catchment in this study has a mean annual precipitation of more than 200 cm. The   Catchment characteristics, including baseline 1980s (1970Catchment characteristics, including baseline 1980s ( -1999 mean annual precipitation (P), mean annual temperature (T), and mean winter (1-Nov to 1-Apr) temperature (T w  headwaters of the Capilano River fill the Capilano Reservoir, which supplies one-third of the water supply for the 2.5 million Metro Vancouver residents (Metro Vancouver 2017).

Methods
Groundwater-surface water (GW-SW) modelling Groundwater discharge during low flow and drought periods is dependent on the amount of snow that accumulates during winter and on the timing and rate of snowmelt and resulting groundwater recharge in the spring and summer (Tague and Grant 2009;Godsey, Kirchner, and Tague 2014;Meixner et al. 2016). Therefore, a comprehensive approach for analyzing climate change impacts on snow drought and streamflow drought requires the application of a distributed, physically-based groundwater-surface water model. With this type of model, parameters are directly related to the physical characteristics of the catchment. Compared to the two other main types of hydrological models, i.e. empirical models and lumped conceptual models, distributed physicallybased models are more appropriate for simulating ungauged catchments (Refsgaard and Knudsen 1996) and for use where significant changes in catchment conditions are expected (Kleme s 1985; Refsgaard and Knudsen 1996), e.g. climate change scenario modelling. Therefore, this study uses the distributed physically-based GW-SW modelling code MIKE SHE -MIKE 11 (Danish Hydraulic Institute (DHI), 2007). MIKE SHE has been used in previous climate change scenario modelling studies, including studies in catchments with seasonal snow cover (Liu et al. 2011;Thompson 2012;Foster and Allen 2015), and has been compared to other modelling codes and shown to adequately model stream discharge (Vansteenkiste et al. 2013;Golmohammadi et al. 2014) and groundwater recharge (Foster and Allen 2015). MIKE SHE is a fully distributed hydrologic modelling code that can simulate actual evapotranspiration (AET), overland flow, one-dimensional (1 D) unsaturated flow, and three-dimensional (3 D) variably saturated groundwater flow. Rivers, lakes, and other channels are simulated by the MIKE 11 model, which is coupled to the MIKE SHE model with river links (h-points). Further details on the comprehensive modelling capabilities of the MIKE SHE software can be found in the user manual (DHI, 2007).
Model boundary conditions were consistent between all models and consisted of: These boundary conditions direct the input (precipitation) onto the model domain and partition it into two outputs: 1) evapotranspiration or 2) surface water flow through the downstream flux boundary. The head-dependent boundary condition applied at the downstream outflow assumes that all water from the catchment exits as streamflow. While groundwater can bypass this discharge point and exit the catchment as deep groundwater flow, the volumetric rate is generally quite low (e.g. Voeckler, Allen, and Alila 2014). Therefore, deep groundwater flow exiting the catchment was excluded from the models. Initial groundwater levels were assigned to coincide with the ground surface and declined to dynamically stable levels during model spin up. The models were run for 150-year periods (1950-2100), using the first 20 years (1950)(1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969) as the spin up period to achieve a dynamically stable state.
The GW-SW models in this study were developed to explore the relationship between snow drought and hydrological drought in idealized systems. Each of these small, ungauged headwater catchment models represent a different ecoregion in British Columbia. While real-world topography, stream networks, etc. were used, all models in this study are generic and were not put through any calibration or validation procedure; therefore, these models provide a basis for comparison (Anderson, Woessner, and Hunt 2015) and not quantitative predictions. Small headwater catchments were chosen due to their hydrological importance (Fassnacht, Webb, and Sanford 2017;Wohl 2017) as well as to limit the computational demand. The modelling code used in this study meets the guidelines laid out by Freeze and Harlan (1969) for adequate physics-based hydrological modelling and the climate change application criteria of Kleme s (1985). Beven (1989), Grayson, Moore, and McMahon (1992), and Voss (2011aVoss ( , 2011b, however, warn against the over parameterization of physically-based models. This study aimed to represent the hydrologic systems in the simplest way possible and used homogenous land cover, soils, and geology for each catchment. Land surface data and overland (OL) flow Actual transpiration and soil evaporation were calculated using the equations from Kristensen and Jensen (1975). Required inputs include leaf area index (LAI), canopy interception, root characteristics, and empirical coefficients. Leaf area index (LAI) for each catchment was estimated using the 10-day interval LAI dataset from Gonsamo and Chen (2014). The required root characteristics include root depth and a root mass distribution parameter Aroot. Aroot was left at the default value (0.25 m À1 ; DHI, 2007). Root depth was set at 400 mm for all catchments based on the findings of Curt, Lucot, and Bouchaud (2001). The canopy interception parameter C int and the empirical coefficients C1, C2, C3 were left at the default values (0.05 mm and 0.3, 0.2, and 20 mm/day, respectively; DHI, 2007), following Voeckler, Allen, and Alila (2014) and Foster and Allen (2015).
Overland flow occurs via two mechanisms: 1) exceedance of the soil infiltration capacity or 2) water table intersects the ground surface and generates return flow. Within MIKE SHE, overland flow is routed by surface topography, with the rate dependent on the diffusive wave approximation of the Saint Venant equation. Topography for all models was assigned using the 0.75 arcsecond Canadian Digital Elevation Model (CDEM; Natural Resources Canada 2017). Land cover data (DataBC 2013) were used to define surface roughness values (Manning's M), which control the resistance to overland flow. Manning's M is the reciprocal of Manning's n; therefore, Manning's n values from Chow (1959) were used to estimate Manning's M values (Table 1).

Unsaturated and saturated zone
The unsaturated zone (UZ) within MIKE SHE is the zone through which the water table rises and falls. Vertical flow in the UZ was modelled using Richards' equation (Richards 1931). Three-dimensional groundwater flow in the saturated zone (SZ) is based on Darcy's equation and is solved implicitly using a finite difference technique. The UZ and SZ are explicitly coupled, and the upper boundary of the SZ is a flux boundary which receives recharge from the SZ. Flux from the unsaturated to the saturated zone varies in time and is computed at the interface of the two zones. Because of the coupling between the two zones, the UZ and SZ must overlap, with the UZ extending to a depth of the lowest possible groundwater head. Depth to groundwater in the catchments is unknown; therefore, to ensure coupling between the UZ and SZ zones, the UZ and SZ zones were assigned the same number of layers, with the same depths, and the same hydraulic properties, extending from ground surface to 200 metres depth.
UZ and SZ layer depths, bulk densities, and vertical and horizontal saturated hydraulic conductivities (K z , K xy ) were assigned based on British Columbia (BC) soil descriptions (Agriculture and Agri-Food Canada 2017). Textural classes were determined from the BC soil descriptions and used to assign parameter values of effective porosity (h s ), residual porosity (h r ), and empirical constants a and n based on values in Carsel and Parrish (1988). Specific yield values were estimated from Morris and Johnson (1967). Parameters for the organic soil layers in the Fort Nelson (Taiga ecoregion) catchment were based on values from Letts et al. (2000). Soil names, bedrock geology, and the associated parameter values are provided in the Appendix.

MIKE 11 stream network
Stream routing was modelled within MIKE 11 and requires four main components: a stream network, cross-sections, boundary conditions, and hydrodynamic parameters. Stream networks and drainage boundaries were obtained from British Columbia Freshwater Atlas (British Columbia 2013). For the stream cross-sections, floodplain topography was extracted from the CDEMs and expert knowledge and hand-digitizing was used to estimate geometry of the main channel. Rating curves for the downstream head-dependent flux boundaries were calculated using Manning's equation, shown in Equation (1).
where Q [m 3 /s] is the discharge leaving the model domain, A [m 2 ] is the cross-sectional area, n is Manning's roughness coefficient, R is the hydraulic radius [m], and S [m/m] is the channel slope. For each catchment, the channel slope was calculated from the CDEM and the hydraulic radius was calculated from the downstream cross-section for a range of possible stream stages spanning low flow to high flow conditions. A value of Q was determined for each stage value using Equation (1), thereby creating the rating curve required for the head-dependent flux boundary. The global bed resistance Manning's roughness coefficient (Manning's n) was set to 0.05 for all catchments based on values in Chow (1959). Conductance values, which control water flow between the stream network and the saturated zone, were estimated from the vertical hydraulic conductivity values of the upper soil horizons (Appendix).

Climate data
Statistically downscaled climate data from three global climate models (GCMs) from Phase 5 of the Coupled Model Intercomparison Project (CMIP5) were used as input to the model. The historical period is represented by the baseline 1980s  data. The models were also run for two future periods, 2050s (2041-2069) and 2080s (2070-2099). Daily climate time series (maximum temperature, minimum temperature, and precipitation) were downloaded from the Pacific Climate Impacts Consortium (PCIC) data portal (PCIC, 2014) covering the period of 1950 to 2100. Mean daily temperature was calculated as the average of the minimum and maximum daily temperature and used as the input for MIKE SHE. A three-member ensemble of GCMs (CNRM-CM5-1, CanESM2-r1, ACCESS1-0-r1) was selected following recommendations in Cannon (2015) with the goal of capturing the widest spread in possible future climate. The GCM data are statistically downscaled with bias-correction/constructed analogues with quantile mapping reordering (BCCAQ). Werner and Cannon (2016) showed that, out of the seven downscaling methods tested, BCCAQ performed best for reproducing hydrologically relevant climate extremes. Two emissions pathways, representative concentration pathways (RCPs) 4.5 and 8.5 were selected for future climate projectionsthese represent different trajectories of anthropogenic radiative forcing, leading to radiative forcing levels of 4.5 and 8.5 W/m 2 by the end of 21st century (van Vuuren et al. 2011). RCP 4.5 represents a medium stabilization scenario, and RCP 8.5 represents a very high baseline emissions scenario (van Vuuren et al. 2011). A comparison of climate data between the baseline 1980s  and two future periods, 2050s (2041-2069) and 2080s (2070-2099) is included in the results.

Evapotranspiration
In addition to mean daily temperature and daily precipitation, MIKE SHE requires estimates of potential evapotranspiration (PET). Potential evapotranspiration (PET) was calculated with the FAO Penman-Monteith method (Allen et al. 1998) using the R package 'sirad' (Bojanowski 2016). Daily solar radiation inputs for the Penman-Monteith method were estimated from daily maximum and minimum temperature using the Hargreaves and Samani (1985) model, following recommendations in Aladenola and Madramootoo (2014). Estimates of daily mean wind speed were unavailable, and a constant wind speed of 5 km/hr was used for all PET calculations, which is within the range of climate normals for the nearby climate stations (Environment and Climate Change Canada 2018). Wind speed exhibits relatively minor impacts on PET (McKenney and Rosenberg 1993; Gong et al. 2006;Tabari and Talaee 2014;C ordova et al. 2015); therefore, the use of a constant wind speed was deemed acceptable.

Snow
Within MIKE SHE, snow accumulation and melt are modelled using a threshold melting temperature, a maximum wet snow storage fraction, and a degreeday coefficient. The threshold melting temperature for all catchments was set to 0 C, and the maximum wet snow storage fraction was set to 0.2, which is in the mid-range of values used in previous studies (Voeckler, Allen, and Alila 2014;Wijesekara et al. 2014;Foster and Allen 2015). A value of 2.74 mm/ degree-day C was used for the degree-day coefficient in all models based on recommendations in United States Department of Agriculture (USDA) National Engineering Handbook (Van Mullem and Garen 2004). The minimum snow storage was set to 0 mm for all catchments.
This snowmelt methodology, referred to as the temperature-index or degree-day method, assumes an empirical relationship between air temperatures and melt rates and has been widely applied due to its simplicity (e.g. Clyde 1931; Corps of Engineers 1956; World Meteorological Organization (WMO), 1986; Van Mullem and Garen 2004). The degree-day method does not, however, account for several factors that are important for snowmelt, including wind speed, humidity, topography (slope, aspect, and shading), cloud cover, and vegetation (Male and Granger 1981;Gray and Landine 1988;Harding and Pomeroy 1996;Pomeroy et al. 1998;Marks et al. 1999;among others). Despite the over-simplification and documented short-comings of this method, temperatureindex methods often perform well at the catchment scale and coarse temporal scales (World Meteorological Organization (WMO), 1986; Sand 1992; Rango and Martinec 1995;Hock 2003). Despite these differences, the relationship between the baseline and future simulations is the samea continuing decline in winter snowpack. While further study using energy balance models should be completed to investigate within-catchment spatial differences in the response of snow drought regimes to climate warming, the degree-day method was deemed sufficient for comparing the average catchment conditions.
Comparison with observed streamflow While no observed streamflow data exist for the four headwater catchments simulated in this study, simulated streamflow for two of the catchments, Capilano and Whiteman, was compared to observed streamflow data from downstream gauging stations as a simple reality check. Simulated streamflow from the Capilano headwater catchment was compared to the observed streamflow from gauge 08GA010 (Capilano River above intake; watershed area 173 km 2 ). Simulated streamflow from the Whiteman headwater catchment was compared to the observed streamflow from gauge 08NM174 (Whiteman Creek above Bouleau Creek, watershed area 114 km 2 ). Comparisons between the Fort Nelson and Blueberry catchments and observed data was not completed due to a lack of appropriate downstream streamflow data. For consistency between the two gauging stations and the model simulations, the 1971-2010 period was used to compare the intra-annual and inter-decadal variability of the simulated and observed streamflow.
To account for the large differences in watershed area (<10 km 2 for simulated catchments versus > 100 km 2 for gauged catchments), the simulated and observed streamflow were transformed from discharge (units of m 3 /day) to runoff (units of mm/day). The time series were then smoothed with a 30-day centered moving average filter to eliminate day-to-day variations. The smoothed runoff was normalized by dividing by the 1971-2010 mean daily runoff, and decadal average hydrographs (1970s, 1980s, 1990s, 2000s) were calculated to analyze the intra-annual and decadal variability of the simulated versus observed streamflow. The two-tailed Spearman's correlation and the volumetric efficiency (VE) were used to quantify similarity between the simulated and observed decadal hydrographs. VE (Criss and Winston 2008) represents the fraction of water delivered at the proper time and is a useful metric for comparing annual hydrographs. VE values range from 0 to 1, where 1 represents a perfect match.

Snow drought
To investigate how different snow drought types impact seasonal low flows, snow droughts were classified using the methodology outlined by Dierauer, Allen, and Whitfield (2019). This classification method uses the mean daily temperature and daily precipitation from the climate models and the simulated snow water equivalent (SWE) data from the MIKE-SHE models. As described in Dierauer, Allen, and Whitfield (2019), this method uses a simplified summer versus winter seasonal classification, which was defined independently for each catchment. For all metrics, the baseline averages were defined using the baseline 1980s period . The winter season was defined using the baseline 25 th percentile of mean daily temperature (T 25 ), with days of the year with T 25 <0 C corresponding to 'winter' and all other days of the year corresponding to 'summer'. Years with below-average peak SWE were identified as snow drought years.
The type of snow drought was classified based on winter precipitation (P w ) and winter thawing degrees (TD w ), where TD w was calculated as the sum of positive mean daily temperatures during the winter season and P w was calculated as the total precipitation during the winter season. Snow drought years with belowaverage P w and below-average TD w are classified as 'dry'; snow drought years with above-average P w are classified as 'warm'; and snow drought years belowaverage P w and above-average TD w are classified as 'warm and dry'.
Following snow drought classification, severity (mean deficit below baseline average), frequency (fraction of years), and risk (frequency x severity) were calculated for the three different snow drought types for each of the climate scenarios. Using the 1980s baseline average peak SWE to define snow drought years results in the identification of many 'minor' snow drought events, where the peak SWE levels are near the baseline average. Warm snow droughts, however, are likely to be minor events and a high threshold was chosen because a lower threshold would likely exclude many temperature-based SWE anomalies (Dierauer, Allen, and Whitfield 2019).
Low flows Low flows are generated by different hydrological processes in winter versus summer, i.e. below freezing temperatures and snow accumulation versus above freezing temperatures and high evapotranspiration rates (Waylen and Woo 1987;Laaha and Bl€ oschl 2006;Burn et al. 2008). Therefore, the same simplified seasonal classification outlined used for the snow drought analysis was applied to the low flow analysis. The magnitude of low flows was quantified using two metrics: 1) mean 15-day minimum runoff (MAM15) and 2) mean 30-day minimum runoff (MAM30). The two metrics were calculated separately for the summer and winter season, indicated by 's' and 'w' subscripts, respectively (Table 2).

Streamflow drought
To evaluate the propagation of snow drought into summer streamflow drought, the frequency and mean severity of summer streamflow droughts following a snow drought were tabulated, considering only years with summer precipitation above the 1980s baseline normal, where 'summer' is defined as earlier in the snow drought classification methods. By doing this, summer streamflow droughts caused by a summer precipitation deficit could be separated from summer streamflow droughts caused by snow droughts, and thus the propagation of snow drought into summer streamflow drought could be highlighted. Summer streamflow drought severity was analysed using the summer low flow metrics, identifying low flows of anomalously low magnitude based on an exceedance threshold (Table 2). A baseline threshold of 80% exceedance frequency was used, i.e. low flow magnitude that was exceeded 80% of the time during the baseline period. Summer streamflow drought years were then defined as years with MAM15s and/or MAM30s values below the threshold.

Water balance analysis
A water balance analysis was completed for each GCM-RCP combination, for a total of six 150-year water balance time series for each headwater catchment. The water balance components of interest for this study include precipitation, snow water equivalent (SWE), runoff, actual evapotranspiration (AET), and groundwater recharge. Groundwater recharge represents water that enters the SZ from the UZ; however, with the MIKE SHE water balance tool, both flux into and out of the SZ are tabulated. Thus, recharge can be positive (outward flux to UZ) or negative (inward flux to SZ). With no change in groundwater storage, the inward fluxes in recharge zones would be balanced by outward fluxes in discharge zones, and total annual recharge would roughly equal 0. Therefore, a detailed saturated zone water balance was extracted for a small upland region of each catchment to analyse groundwater recharge in a recharge zone only. Runoff for each catchment was calculated as the volume of streamflow leaving the downstream head-dependent flux boundary, converted from m 3 /s to mm/day by dividing by catchment area.
Water balance results for individual GCMs are not shown but rather lumped by RCP for simplicity. Results are summarized for a baseline period (1980s) and two future periods (2050s and 2080s) and are presented as both the absolute change (futurebaseline) and relative change ([futurebaseline]/baseline) for the three-member GCM ensemble (Figure 3 and Appendix). The average annual numerical water balance errors, which stem from the uncertainty in algorithms used for the finite difference modelling methods, were less than 3% for all models. However, due to transient conditions and changes in subsurface storage, mean annual AET plus mean annual runoff is not equal to mean annual precipitation.

Simulated versus observed streamflow
Modelled streamflow from the uncalibrated Capilano and Whiteman headwater catchment exhibits seasonal patterns and inter-decadal variability that resembles observed patterns in streamflow at the downstream gauging stations (Figure 2). Correlation between the simulated and observed streamflow is strong (Spearman's rho ¼ 0.61-0.65) and similarity in the seasonal timing of observed versus simulated streamflow is strong for the Capilano catchment (VE ¼ 0.66) and moderate in the Whiteman catchment (VE ¼ 0.48). For both catchments, the modelled streamflow exhibits an earlier freshet peak compared to the observed data. This difference in freshet timing is likely due to scale effects, as the observed data are from catchments that are much larger than the simulated headwater catchments (>100 km 2 versus <10 km 2 ). The gauged catchments include heterogenous land use/landcover, soil properties, geology, etc., along with regions of higher altitude and latitude compared to the small headwater catchment which were modelled with homogenous catchment properties. Thus, differences in freshet timing are expected, and this comparison with observed data demonstrates the ability of the GW-SW models to capture the intra-annual and decadal streamflow variability.

Climate change impacts on the annual and intraannual water balance
The statistically downscaled climate change projections show increases in temperature and precipitation for all four catchments (Figure 3 and Appendix). The relative increase in precipitation is highest in the two northern catchments (Fort Nelson and Blueberry), which have a cold, dry climate. The absolute increase (mm/year) in precipitation, however, is highest in the warmest, wettest catchment -Capilano. The seasonal distribution of precipitation does not change substantially under RCPs 4.5 and 8.5; however, in general, the fall (Sep to Nov) and spring (Mar to May) seasons exhibit the largest relative increases in precipitation (Figure 3), and summer (Jun to Aug) and winter (Dec to Feb) exhibit the smallest. Compared to the projected changes in precipitation, changes in temperature are more similar among all catchments (Appendix). Projected temperature increases are highest in winter and lowest in fall ( Figure 4). As expected, increases in the mean annual temperature are greatest for the high emissions scenario (RCP 8.5), and by the 2080s, mean annual temperature is projected to be 5.6 to 6.1 C higher compared to the baseline 1980s period (Appendix). Annual runoff is projected to increase for all four catchments ( Figure 3). Relative (%) increases in runoff are largest for the coldest, driest catchment (Fort Nelson), while absolute (mm/year) increases in runoff are largest for the warmest, wettest catchment (Capilano; Figure 3 and Appendix). In addition to increases in annual runoff, the within-year distribution of runoff changes substantially. The spring freshet starts earlier for all catchments and decreases in magnitude for all but Fort Nelson (Figure 4), which is the northernmost and coldest catchment. In the warmest catchment (Capilano), the spring freshet disappears completely for both future time periods under both RCPs. In all catchments, the slope of the spring freshet rising limb decreases, indicating a longer spring melt season with slower snowmelt (Figure 4). These changes (declined spring freshet peak and a longer melt season) are consistent with the findings of previous observation-based studies (Hamlet and Lettenmaier 2007;Rood et al. 2008).
In general, changes in the timing and magnitude of the spring freshet are directly related to changes in the length of the snow-covered period and the magnitude and timing of peak SWE (Figure 4). While the increased annual precipitation leads to increased peak SWE for the coldest catchment (Fort Nelson), the impact of increased temperatures outweighs the impacts of increased precipitation in the remaining catchments, leading to no significant change in peak , and (f) annual groundwater recharge. Shading indicates a significant (p < 0.05) increase or decrease relative to the baseline period, as assessed with the two-sided Mann-Whitney U test. Table A2 in the Appendix provides the corresponding mean annual values along with the absolute and relative change. SWE (Blueberry) or significant decreases in peak SWE (Whiteman and Capilano; Figure 3). The largest absolute and relative decreases in peak SWE occur in the warmest catchment, Capilano, which has a > 90% decrease in peak SWE for the 2080s under RCP 8.5. In addition to changes in the magnitude of peak SWE, the average day of peak SWE and melt-out occur earlier, resulting in a shorter snow-covered period for all catchments for both future time periods under both RCPs (Figure 4). Additionally, snowmelt is slower for all catchments for both future periods, as illustrated by the shallower slope of the falling limbs ( Figure 4). The earlier and slower snowmelt is consistent with the Musselman et al. (2017) study, which showed that snow melts more slowly in a warmer world due to an increase in winter and spring melt and longer snow-free periods during times of high energy (i.e. summer).
As expected, the changes in temperature and precipitation, and associated changes in snow accumulation and melt, lead to significant changes in both the total annual AET (Figure 3) and the intra-annual patterns in AET (Figure 4). AET increases significantly for all catchments except Whiteman, located in the Okanagan Valley, which exhibits no significant change in annual AET (Figure 3). Seasonally, AET increases most in late-winter and spring (Feb to May) and decreases (Whiteman and Capilano) or exhibits no substantial change (Fort Nelson and Blueberry) during summer (results not shown). Decreased AET during the summer can be primarily attributed to the shift toward earlier snowmelt, which decreases summer water availability. This negative feedback between snowmelt timing and evapotranspiration has been discussed by Barnett, Adam, and Lettenmaier (2005) and documented by previous climate change modelling studies (e.g. Shrestha et al. 2012). Shifts in vegetation patterns will likely influence catchment response to climate change (Alo and Wang 2008;Teutschbein et al. 2018); however, it is difficult to project and constrain possible vegetation shifts, and vegetation change was not included in the modelling efforts. Changes in wind speed, either due to climate change or vegetation change, would also impact AET, an effect which was not considered in this study.
Within MIKE SHE, water that reaches the saturated zone (i.e. groundwater recharge) may then exit the saturated zone via evapotranspiration, baseflow to the river, or surface return. Groundwater recharge may be higher or lower than runoff, depending on catchment's physical properties (e.g. soils, geology, vegetation, ground slope) which control the evapotranspiration dynamics and the magnitude of overland flow. In the Fort Nelson catchment, which has high porosity organic soils (Appendix) and shallow slope (Table 1), groundwater recharge is much higher than runoff (Appendix), and a large proportion of the water that reaches the saturated zone then leaves the system through evapotranspiration. In the Capilano catchment, which has lower porosity soils (Appendix) and steep slopes (Table 1), groundwater recharge is lower than runoff (Appendix) and a substantial portion of runoff is generated from overland flow.
At the annual time scale, recharge increases significantly for all but the warmest catchment (Capilano; Figure 3). Intra-annually, the patterns in groundwater recharge are primarily affected by changes in the onset of snow accumulation and melt. In all catchments, the spring recharge peak starts earlier in the year and decreases in magnitude (Figure 4), resulting in higher winter groundwater storage, an earlier start to the spring/summer groundwater recession period, and thus decreased summer groundwater storage (Figure 4). Increased winter-season recharge for regions with seasonal snow cover is consistent with the results of previous climate change modelling studies (e.g. Eckhardt and Ulbrich 2003;Jyrkama and Sykes 2007;Kovalevskii 2007). A shift toward more rain and less snow in combination with slower snowmelt would suggest an overall decrease in groundwater recharge (Earman et al. 2006;Barnhart et al. 2016). While it is difficult to separate the effects of increased temperatures from the effects of increased precipitation, the results of this study show an increasing ratio of recharge to precipitation (R:P ratio) for the Fort Nelson catchment, relatively constant R:P ratios for the Blueberry and Whiteman catchments, and a decreasing R:P ratio for the Capilano catchment (Appendix). The different responses of the R:P ratio (increase, no change, decrease) appear related to the catchment's starting point (in terms of temperature), as the coldest catchment exhibits an increase in the R:P ratio and the warmest catchment exhibits a decrease in the R:P ratio.

Snow drought
In response to the increased precipitation and temperature (Figure 3), the snow drought regime changes substantially for all catchments. Warm snow droughts increase in frequency and dry snow droughts decrease in frequency ( Figure 5). Additionally, warm, and warm and dry, snow drought severity increases for the two warmest catchments, Whiteman and Capilano ( Figure 6). In general, dry snow droughts transition to warm and dry snow droughts, and, by the 2080s, the frequency of dry snow drought drops to 0 for all catchments ( Figure 5). In terms of temperature, the magnitude of change in the snow drought regime is related to the catchment's starting point, with the warmest catchment (Capilano) exhibiting the largest increase in the frequency and severity of snow drought and the coldest catchment (Fort Nelson) exhibiting no substantial increase in the frequency or severity of snow drought (Figures 5 and 6).
The increased frequency and severity of snow drought necessarily leads to increased snow drought risk, and overall, the changes in snow drought risk (Table 3) mirror the changes in snow drought severity ( Figure 6). In general, snow drought regimes in all catchments shift toward more frequent, higher severity warm, and warm and dry, snow droughts, and less frequent, lower severity dry snow droughts. As documented by Dierauer, Allen, and Whitfield (2019), the response of warm snow drought risk to increased winter temperature is non-linear. A 2 C increase in the mean winter (1-Nov to 1-Apr) temperature corresponds to a substantially larger increase in warm snow drought risk for the Capilano catchment as compared to the Fort Nelson catchment. The two warmest catchments, Whiteman and Capilano, exhibit the largest increases in total snow drought risk. Due to the transition of dry snow droughts to warm and dry snow droughts, dry snow drought risk decreases for all catchments for both future time periods under both RCPs. The coldest catchment, Fort Nelson, exhibits a slight decrease in total snow drought risk.

Low flows
As the snow drought regime shifts toward more frequent, higher severity temperature-related (i.e. warm, and warm and dry) snow droughts, the streamflow regime shifts toward more severe summer low flow periods and less severe winter low flow periods (Figure 7). Low flows are an annually recurring feature of the natural flow regime (Smakhtin 2001); however, more severe low flow periods are equivalent to streamflow droughts. Thus, a shift toward more severe summer low flows represents an increase in summer streamflow drought severity, and a shift toward less severe winter low flows represents a decrease in winter streamflow drought severity.
The impact of snow drought on summer and winter low flows is dependent on the snow drought type. Consistent with findings of Dierauer, Whitfield, and Allen (2018), compared to years with above-average peak SWE, years with warm, and warm and dry, snow droughts lead to more severe summer low flows and less severe winter low flows (Figure 8). The impacts of snow drought on summer and winter low flows are significant in the two warmer catchments (Whiteman and Capilano; Figure 8). In the colder catchments (Fort Nelson and Blueberry), however, the impact of snow drought on low flows are more varied. Warm and dry snow droughts, for example, have no significant impact on summer low flow severity in the colder catchments, possibly due to the counteracting effect of increased groundwater recharge and groundwater storage levels during the winter season (Figure 4) which provides more baseflow to sustain higher low flows during the summer. Alternatively, differences in low flow response to the snow drought types may be related to catchment properties. Further research is needed to address the role of catchment control on streamflow response to snow drought, as it is beyond the scope of this study.

Summer streamflow droughts
In the context of climate warming and considering the relationship between snow drought and low flows shown in Figure 8, summer streamflow drought regimes are likely to shift toward more frequent, higher severity snow-drought related events. Using a threshold-based approach to define summer streamflow drought years shows that, in the absence of summer precipitation deficit, snow droughts propagate into summer streamflow droughts more frequently in the future time periods as compared to the baseline 1980s ( Figure 9). Thus, warm, and warm and dry, snow droughts not only become more frequent and severe in the future (Figures 5 and 6) but are also more likely to result in summer streamflow drought conditions. Dry snow droughts, on the other hand, become less frequent in the future ( Figure 5) and, in the absence of summer precipitation deficit, are unlikely to lead to summer streamflow droughts (Figure 9).
To determine if changes in the summer streamflow drought regime are related to baseline climate conditions, the streamflow drought events identified in Figure 9 were lumped together. The frequency of the events was then plotted against the mean winter (1-Nov to 1-Apr) temperature ( Figure 10). This analysis shows that the magnitude of change in the frequency of these streamflow drought events varies between catchments and shows a strong relationship with the baseline winter air temperature. The Fort Nelson catchment, which has a baseline winter air temperature far below zero, exhibits minimal increase (<10%) in the occurrence of snow-drought-related summer streamflow droughts, while the Capilano catchment, which has a baseline winter air temperature near zero, exhibits a large increase (> 40%) in the Figure 6. Mean severity (negative SWE anomaly represented as the fraction below baseline average) of dry (D), warm (W), and warm and dry (W&D) snow droughts for the baseline 1980s versus 2050s (2040-2069) and 2080s (2070-2099) for representative concentration pathway (RCP) 4.5 and RCP 8.5. Note: Dry snow droughts transition to warm and dry snow droughts and therefore have no mean severity plotted in some future time periods due to non-occurrence. Table 3. Risk (severity x frequency) for dry (D), warm (W), and warm and dry (W&D) snow droughts. Baseline 1980s (1970Baseline 1980s ( -1999Baseline 1980s ( ) versus 2050s (2040Baseline 1980s ( -2069Baseline 1980s ( ) and 2080s (2070Baseline 1980s ( -2099 for representative concentration pathways (RCPs) 4.5 and 8.5. occurrence of these summer streamflow droughts ( Figure 10).

Relation to previous work
The results of this study reinforce the findings of previous research on the climate change impacts on snow and its hydrological impacts (Cooper, Nolin, and Safeeq 2016;Sproles, Roth, and Nolin 2017;Hatchett and McEvoy 2018;Marshall et al. 2019). Snow droughts have direct impacts on summer low flows (Figure 8), and temperature-related (i.e. warm, and warm and dry) snow droughts become more frequent and severe and thus more likely to result in summer streamflow drought conditions (Figure 9). The shift to lower severity winter streamflow droughts and higher severity summer streamflow droughts is consistent with the results of previous hydrologic modelling studies (Feyen and Dankers 2009;Wanders and Van Lanen 2015) and with the general hypothesis that streamflow droughts with different causative factors will respond differently to climate change (Van Loon and Van Lanen 2012;Van Loon et al. 2015). In general, increased summer streamflow drought severity and decreased winter streamflow drought severity are consistent with an overall shift in the intra-annual distribution of runoff, an impact of climate warming on snowmelt hydrology which has been documented by many previous studies (Leith and Whitfield 1998;Whitfield and Cannon 2000;Adam, Hamlet, and Lettenmaier 2009;D ery et al. 2009;Pederson et al. 2011;among others). Consistent with the results of Dierauer, Allen, and Whitfield (2019), the response of snow drought risk to climate warming is non-linear, and the magnitude of change in the snow drought and summer streamflow drought regime is related to the catchment's baseline mean winter (1-Nov to 1-Apr) temperature (T w ) ( Figure 10). Because of the nonlinear relationship between temperature and frequency of snow drought propagation to summer streamflow drought, a þ 2 C change in the mean winter temperature has a larger impact on the summer streamflow drought regime in catchments with winter temperatures near  (1970-1999) versus 2050s (2040-2069) and 2080s (2070-2099) for representative concentration pathway (RCP) 8.5. Shading indicates a significant (p < 0.05) increase or decrease relative to the baseline period, as assessed with the two-sided Mann-Whitney U test. Abbreviations are as in Table 2. Figure 8. Snow drought impacts on low flows by snow drought type, including years without snow drought (None) and years with dry (D), warm (W), and warm and dry (W&D) snow droughts. Shading indicates the values are significantly (p < 0.05) higher or lower relative to years without snow drought, as assessed with the two-sided Mann-Whitney U test. Abbreviations are as in Table 2. Figure 9. Frequency of snow drought propagation into summer streamflow drought for RCP 8.5, in the absence of summer precipitation deficit, by snow drought type: dry (D), warm (W), warm and dry (W&D). Abbreviations are as in Table 2. Note: Upper limit of y-axis is 0.5 (50%). Results for RCP 4.5 not shown.
zero (e.g. the Capilano catchment) compared to catchments with winter temperatures far below zero (e.g. the Fort Nelson catchment).

Implications for ecosystem health
The shift toward more frequent and more severe temperature-related snow droughts (i.e. warm, and warm and dry) and more severe summer low flow periods will have wide-ranging effects on terrestrial and aquatic ecosystems. In the terrestrial realm, earlier snow disappearance will lead to less water for mountain ecosystems (Harpold 2016), lower carbon uptake (Hu et al. 2010;Winchell et al. 2016), more wildfires (Westerling et al. 2006), and more tree death (Bales et al. 2018). In the aquatic realm, climate warming coupled with the shift toward longer more severe summer low flows will have compound impacts on stream temperature, which increases with air temperature and decreases with streamflow rates (Hockey, Owens, and Tapper 1982;Webb, Clack, and Walling 2003). Summer low flow periods are critical for aquatic ecosystem health (Fleming et al. 2007;Moore, Nelitz, and Parkinson 2013), and higher stream temperatures negatively impact species distributions and decrease growth rates (Beschta et al. 1987, Eaton andScheller 1996). Further, summer low flows that are lower than normal, i.e. drought conditions, reduce habitat availability (Lake 2003), increase pollutant concentrations (Mosley 2015), and lower the oxygen available to aquatic organisms (Sprague 2005;van Vliet and Zwolsman 2008;Ylla et al. 2010).
Shifts in snow drought, low flows, and streamflow drought regimes will also have widespread implications for surface water supply security. Increased frequency of warm snow droughts will likely lead to an increased frequency of mid-winter melt events (Hatchett and McEvoy 2018), which will create challenges for reservoir management. Winter melt events should be of low intensity (e.g. Musselman et al. 2017); however, climate change may also result in increased rain on snow events and thus high intensity flows, i.e. floods (Yan et al. 2018). Reservoirs may require higher flood control due to the increased winter flows and increased risk to rain on snow events, while simultaneously requiring more storage capacity to counter decreasing summer flows. As summer low flow periods become more severe and snow-drought related summer streamflow droughts become more frequent, summer water scarcity may increase. The most severe water scarcity will likely occur due to the coincidence of warm and dry conditions (AghaKouchak et al. 2014) and layered impacts from different drought types (Van Loon et al. 2015).

Study limitations and implications for future water management
The GW-SW models in this study are generic, and, therefore, represent interpretive tools (Anderson, Figure 10. Frequency of snow drought propagation into summer streamflow drought versus the mean winter (1-Nov to 1-Apr) temperature (Tw) for the historical and future climate scenarios. Linear trend lines plotted by watershed. Locally weighted smoothing (LOESS) lines (solid black) plotted for all data points to highlight the non-linear relationship between drought frequency and Tw.
Woessner, and Hunt 2015) rather than predictive tools. Like all GW-SW models, they are simplified numerical representations of natural flow systems and cannot duplicate the natural flow system exactly. The models are physically based, however, and the simulated streamflow regime is similar to that observed in downstream watersheds (Figure 2). Further, the GW-SW models simulated the response of low flows and streamflow droughts to snow drought in a manner that is consistent with previous observation-based Safeeq et al. 2013;Jenicek et al. 2016;Kormos et al. 2016;Cooper et al. 2018) and model-based Tague and Grant 2009;Van Loon et al. 2015) studies. Changing input parameter values and incorporating freeze/thaw soil dynamics into the GW-SW models could either decrease or increase sensitivity of snow drought and low flows to climate warming. Selecting catchments across different ecoregions allows broader patterns to emerge, and the catchment's baseline climate conditions likely outweigh parameter uncertainty and model limitations in terms of the relative sensitivity to warming between these four catchments, as suggested by the dominance of climate over catchment control shown by Beck, De Roo, and van Dijk (2015). Thus, the results highlight how a catchment's starting point (i.e. the baseline climate) has a large impact on how sensitive it is to climate warming. The results of this study can provide general insights into future water management challenges and could also be used as a base for identifying areas of interest and designing comparative snow drought and streamflow drought modelling studies. In general terms, the results can be applied to regional water management challenges, as in the following paragraphs.
In NEBC, where the Fort Nelson and Blueberry catchments are located, shifting snow and streamflow drought regimes will likely lead to decreased freshwater security. The unconventional oil and gas industry in NEBC expanded rapidly from 2005-2015, and short-term requirements for large quantities of water for hydraulic fracturing operations (Rivard et al. 2014) put high water demands on local watersheds. In 2015, 45% of the 7.74 million m 3 of water used for hydraulic fracturing in NEBC was sourced from surface water (British Columbia Oil and Gas Commission 2016), with the highest water demand occurring in the warmer months from May to September. Without significant commitment on the part of industry to re-use and recycle water for hydraulic fracturing, industrial water demand is likely to increase substantiallywith the possibility of a more than 350% increase by 2030 compared to 2015 levels under a high development scenario (Kniewasser and Horne 2015). Industrial freshwater abstractions are suspended during streamflow drought conditions, and the British Columbia Oil and Gas Commission has issued water use suspensions four times in the last eight years (https://www.bcogc.ca/directives). As summer low flows decrease, water use suspensions are likely to become more frequent, and balancing increasing demand with decreasing security will be a significant challenge for the region in the future.
The Whiteman catchment is located in the Okanagan Valley, where surface water sources supply 67% of the annual water demand (Summit Environmental Consultants Inc 2010). Most streams in the Okanagan are fully allocated, with no further allocations possible (Cohen and Kulkarni 2001). The greatest proportion water is used for agriculture, and irrigation, which accounts for 75% of regional consumptive water use, is expected to increase considerably with continued climate warming (Neilsen et al. 2006). Additionally, average per person water use is high (Summit Environmental Consultants Inc 2010) and population is expected to grow at a rate of 0.2 to 0.7% per year (BCstats 2017). Population growth in the Metro Vancouver region, where the Capilano catchment is located, is expected to be even higher at 0.9 to 1.4% per year (BCstats 2017). Significant opportunities exist for demand-side reductions in water use for the Okanagan (Neale, Carmichael, and Cohen 2007;DHI Water and Environment 2010) and Metro Vancouver (Metro Vancouver 2011) regions. Water shortages have already occurred both regions (Okanagan Water Stewardship Council 2008;Metro Vancouver 2015), and, considering the results presented in this study and others (DHI Water and Environment 2010;Harma, Johnson, and Cohen 2012), summer water shortages are likely to be more common in the future.
The results presented provide a basis for comparison of snow droughts and their impact on groundwater recharge from snowmelt that provides insight that can be used by water managers to anticipate the types of changes that will occur in British Columbia that will need to be accommodated in future planning. The results presented here should not be used in water allocation since the absolute magnitude and timing of the changes should be addressed with a full energy balance model; however, the results show that, despite increasing precipitation, climate warming can be expected to lead to increased frequency and severity of temperature related snow droughts across four ecoregions in British Columbia.

Conclusion
Climate change impacts on snow drought, low flows, and summer streamflow drought were investigated using generic coupled GW-SW models for four headwater catchments in British Columbia. Results show that increased precipitation and temperature lead to decreased risk of dry snow droughts and increased risk of warm, and warm and dry, snow droughts. Climate warming and shifting snow drought regimes result in decreased summer runoff, decreased summer groundwater storage, and longer and more severe summer low flow periods. Snow droughts have direct impacts on summer low flows, and temperaturerelated snow droughts not only become more frequent and severe in the future but are also more likely to result in summer streamflow drought conditions.
The response of snow hydrology to climate warming is non-linear, and catchments with winter temperatures near 0 C exhibit substantially larger impacts from þ2 C of warming compared to catchments with winter temperatures far below 0 C. The shift toward more frequent and more severe temperature-related snow droughts will decrease water availability during the summer for agricultural and industrial usespotentially leading to decreased freshwater supply security, and the increased frequency of warm snow droughts will likely lead to an increased frequency of mid-winter melt events that will affect reservoir management. Changes in the low flow regimes driven by a warming climate will affect the ecology of river systems.