Assessment of Future Risks of Seasonal Municipal Water Shortages Across North America

While anthropogenic climate change poses a risk to freshwater resources across the globe through increases in evapotranspiration and temperature, it is essential to quantify the risks at local scales in response to projected trends in both freshwater supply and demand. In this study, we use empirical modeling to estimate the risks of municipal water shortages across North America by assessing the effects of climate change on streamflow and urban water demand. In addition, we aim to quantify uncertainties in both supply and demand predictions. Using streamflow data from both the US and Canada, we first cluster 4,290 streamflow gauges based on hydrograph similarity and geographical location. We develop a set of multiple linear regression (MLR) models, as a simplified analog to a distributed hydrological model, with minimum input data requirements. These MLR models are calibrated to simulate streamflow for the 1993–2012 period using the ERA5 climate reanalysis data. The models are then used to predict streamflow for the 2080–2099 period in response to two climate scenarios (RCP4.5 and RCP8.5) from five global climate models. Another set of MLR models are constructed to project seasonal changes in municipal water consumption for the clustered domains. The models are calibrated against collected data on urban water use from 47 cities across the study region. For both streamflow and water use, we quantified uncertainties in our predictions using stochastic weather generators and Monte Carlo methods. Our study shows the strong predictive power of the MLR models for simulating both streamflow regimes (Kling-Gupta efficiency >0.5) and urban water use (correlation coefficient ≈0.7) in most regions. Under the RCP4.5 (RCP8.5) emissions scenario, the West Coast, the Southwest, and the Deep South (South-Central US and the Deep South) have the highest risk of municipal water shortages. Across the whole domain, the risk is the highest in the summer season when demand is high. We find that the uncertainty in projected changes to the water demand is substantially lower than the uncertainty in the projected changes to the supply. Regions with the highest risk of water shortages should begin to investigate mitigation and adaptation strategies.


INTRODUCTION
Based on current climate projections, global mean temperature is likely to increase by at least 2°C, relative to the 1850s, by the end of the century (Forster et al., 2020). This warming will lead to unprecedented consequences for life on this planet including increased wildfire risk (Goss et al., 2020), decreased snowpack (Ashfaq et al., 2013), and decreased freshwater biodiversity (Reid et al., 2019). While every location on Earth is expected to experience increases in near-surface air temperature, ranging from a 1°C increase over oceans to over 5°C in higher latitude areas, the projected changes in precipitation are less certain and less spatially uniform (Hoegh-Guldberg et al., 2018). The projections of hydrologic drought and floods are consequently uncertain, while at the same time, highly costly in terms of their impact on property damage, food shortage, loss of jobs, and loss of lives (Howitt et al., 2015;Achakulwisut et al., 2018;Tellman et al., 2020). More accurate predictions of these extreme events can lead to better mitigation and adaptation procedures, such as green infrastructure, placing restrictions on water consumption, moving water in above-ground reservoirs to below-ground aquifers, or investing in technologies that could improve water use efficiency (Tanaka et al., 2006;Mei et al., 2018;Yang and Liu, 2020).
As anthropogenic climate change progresses and populations across North America continue to grow, freshwater resources on the continent may experience more strain. The United States (US) and Canada use the most water per capita (>1,000 cubic meters per person per year) compared to other G8 countries (100-900 cubic meters per person per year) (Ritchie and Roser, 2017). Despite the use of less than 20% of their available water resources every year, freshwater resources and potential vulnerability to water shortages are non-uniformly distributed throughout the continent (Rosegrant and Cai, 2002). Such shortages will act as a vulnerability multiplier leading to socioeconomic and physical health deterioration in groups such as migrants, poor families in urban areas, and farmers (Sugg et al., 2020). The Canadian Prairies are known to have been vulnerable to historical hydrologic droughts, although few assessments of drought risk across Canada have been made to date (MacDonald et al., 2008;Bonsal et al., 2011). In the US, the regions most severely hit by recent droughts include the Central US (Basara et al., 2019), California (Howitt et al., 2015), as well as the relatively wet areas in the US South . Further droughts in these areas as well as more widespread droughts throughout North America may come as a result of declining summer precipitation in the latter half of the 21st century, warming Pacific and North Atlantic Oceans, and escalating climate variability (Rosegrant and Cai, 2002;MacDonald et al., 2008;Schwalm et al., 2012). Droughts, however, are not the only climate pattern of concern. Snow is crucial in the western United States for sustaining water demand, thus decreasing snowpack as a result of increasing temperatures threatens water sustainability, though uncertainties remain large (Mankin et al., 2015).
Risk of freshwater shortage is measured in a variety of ways (Rosegrant and Cai, 2002;Foti et al., 2012;Dickson and Dzombak, 2019). For example, Foti et al. (2012) defined shortage risk as the difference between water supply and demand, or more specifically, as a probability for crossing the critical threshold when water demand exceeds supply. Different types of process-based or mechanistic models (e.g., Chien and Knouft, 2013;Mahat et al., 2017) and statistical or purely datadriven models (e.g., Barbarossa et al., 2017;Mendoza et al., 2017;Brunner et al., 2020) have been used to project changes in water supply and water demand globally and in North America. Mechanistic models are efficient tools for prediction of water supply and demand in regions where a wide range of data on climatic, landscape, and socioeconomic and demographic attributes are available. Despite great progress being made on developing advanced mechanistic models (Chen et al., 2017), the application of these models in ungauged watersheds, which cover more than 90% of lands globally and across North America, is unpractical due to poor data availability required for building such models (Blöschl et al., 2013). Statistical models, on the other hand, have become known as simple and fast tools for providing general insights on the estimation and forecast of water supply and demand with minimum data requirements. For example, simple regression techniques have been used to relate physical and climate characteristics to hydrological signatures such as flow duration curves and low flow statistics (Jehn et al., 2020). While statistical models are limited in their interpretation of causality among variables and rely on assumptions that cannot be adequately tested, these simple models are widely used to provide a general picture of ungauged watersheds' responses to climate change (Besaw et al., 2010;Razavi and Coulibaly, 2013;Saadi et al., 2019), and informing water resources planning and management. Further, regression models can capture implicit relationships between runoff and explanatory variables for which there is no theoretical explanation due to the coevolution of climate, geology, and topography (Blöschl et al., 2013).
Most assessments of water supply and demand management focused on local areas across North America and only study either supply or demand with regression or physically based models (Balling and Gober, 2007;Breyer and Chang, 2014;Shamir et al., 2015;Byun et al., 2019). Foti et al. (2012), however, provided a comprehensive risk assessment of water shortages for the contiguous US using regression and physically based models. They projected annual water supply and consumption in response to changes in temperature and precipitation from three Global Climate Models (GCMs). Areas such as Central California, Utah, Arizona, and Central US were found to be most vulnerable to water shortages. A more recent assessment of future changes in streamflow regimes across the US found that the changes are expected to be most pronounced for currently melt-dominated regimes in the Rocky Mountains (Brunner et al., 2020). While five GCMs were used in the study, very little agreement was found among the GCM projections, pointing to a relatively large uncertainty caused by the choice of GCMs. In terms of changes to water consumption across the US, negligible changes were found when the projections were based solely on population growth, income growth, and changes in water-use efficiency, i.e., without accounting for climate change (Rosegrant and Cai, 2002;Foti et al., 2012). With climate change remaining a factor, however, water consumption was projected to substantially increase, mainly due to expected increases in agricultural and landscape irrigation in response to rising potential evapotranspiration (Foti et al., 2012;Brown et al., 2013).
As the population and suburb-centered development continue to grow in North America, many municipalities in arid and semiarid regions are acquiring water rights from agriculture in anticipation of an uptick in municipal water demand (MacDonald, 2010;Sabo et al., 2010;Maas et al., 2017). The percent of water demand used for domestic purposes is set to increase, thus the transfer of water resources from agricultural to municipal use represents an additional stressor for rural sustainability and food security (Rosegrant and Cai, 2002). Thus, changes in municipal water use in response to climatic and non-climatic factors are expected to have implications for regional water consumption as a whole. In the US, residential water use accounts for ∼60% of municipal water use (Dieter et al., 2018), with 22-65% of the residential water use coming from outdoor water use (DeOreo et al., 2016), which is shown to be highly sensitive to changes in weather and climate (Gober et al., 2016). This high sensitivity of outdoor water use to climatic drivers has triggered the use of regression models to predict the water consumption based on a combination of weather predictors, such as air temperature, wind speed, precipitation, and evapotranspiration (Wong et al., 2010;Adamowski et al., 2012;Bakker et al., 2014;Chang et al., 2014). The models have revealed a high spatial heterogeneity in the sensitivity of municipal water consumption to changes in these predictors (Opalinski et al., 2020). Maximum daily temperature was found to be the predictor with the most explanatory power across the region, with stronger predictive power in dry regions (Opalinski et al., 2020). In the same study, dry regions were found to have greater seasonality in municipal water use perhaps due to increased seasonal climate variability, with the peak use in summer when there is an increased need for irrigation of the excising green urban landscapes.
Knowledge on the risks of municipal water shortages is essential to inform science-based strategies and decisionmaking tools for water security (Byun et al., 2019;Dilling et al., 2019;Özerol et al., 2020). Current research has focused on forecasting the magnitude and seasonality of water supply (e.g., Chien and Knouft, 2013;Demaria et al., 2016;Mahat et al., 2017;Byun et al., 2019;Brunner et al., 2020) as well as magnitude of water demand (e.g., Ruth et al., 2007;DeOreo et al., 2016;Maas et al., 2017;Opalinski et al., 2020;Rasifaghihi et al., 2020). However, the magnitude and seasonality of water shortages, as a combined effect of water supply and water demand, are still relatively poorly characterized and forecasted in most parts of the globe including North America. Additionally, the forecast of water shortages, similar to other types of forecasts in environmental science, is prone to a wide range of uncertainties, which are generally challenging and difficult to quantify. Such uncertainties are mainly stemmed from uncertainties in time-varying climatic attributes (Katz, 2002), uncertainties associated with more static physical attributes of watersheds such as geology, land cover, and soil (Nilsson et al., 2007;Addor et al., 2017), and uncertainties corresponding to future population growth and water use (Yang et al., 2016;Hart and Halden, 2019;Keilman, 2020). Furthermore, the more complex the modelling framework is, often relying on heavily calibrated mechanistic models of water supply and demand, the more difficult it is to track and quantify its sources of uncertainties (e.g., Foti et al., 2012). On the other hand, simple statistical models with low degrees of freedom present more readily available tools for quantifying these uncertainties. Characterizing and quantifying the risks of future seasonal water shortages, and their associated uncertainties, with the use of simple statistical models were the main motivations of this study.
This study has two main goals: 1) to investigate the use of simple statistical models, based on widely available climatic data, in simulating present and future (end of the 21st century) water shortages under a changing climate across North America (US and Canada); and 2) to assess uncertainties in the projections of water shortages originating from the uncertainties in climate scenarios from the ensemble of GCMs. To address goal 1), multivariate regression models for both supply and demand are developed and calibrated with daily time series of streamflow from river stations and daily time series of municipal water consumption from cities across the region. The models are forced by state-of-the-art climate reanalysis data representing the present climate and weather attributes (e.g., temperature, precipitation, windspeed, evapotranspiration, and rainfall/snowfall characteristics), while the future climate is represented by two emission scenarios from the ensemble of five GCMs. To address goal 2), we use stochastic weather generators as it is one of the recommended methods for estimating uncertainties in climate change projections (Katz, 2002) and in hydrological forecasting (e.g., Caron et al., 2008;Breinl, 2016). By combining the projections of water supply and demand, we aim to provide a spatial pattern of seasonal municipal water shortage risks, together with estimated credibility in these risks across the whole region. 1993-2012. Hourly time series of total precipitation, wind speed, air temperature, dewpoint temperature, and incoming solar radiation are extracted. Dewpoint temperature and air temperature are used to estimate relative humidity with the August-Roche-Magnus approximation (Lawrence, 2005). Cubic inverse distance weighting interpolation is performed to obtain a single time series of each weather variable at each streamflow station. Interpolated precipitation totals that are below 1 mm are set to 0 mm following (Islam and Cartwright, 2020) and (Tian et al., 2009). Evapotranspiration is calculated from the Penman-Monteith equations, as outlined by (Walter et al., 2000) using maximum and minimum air temperature, maximum and minimum relative humidity, wind speed, and incoming solar radiation at the surface. Hourly time series of total precipitation, wind speed, and air temperature at the surface are also extracted at the grid square closest to each city with water consumption data.

Future Climate
Future climate scenarios are obtained from five GCMs from Coupled Model Intercomparison Project phase five [CMIP5, Table 1, (Taylor et al., 2012)], henceforward abbreviated as: GFDL-CM3, HadGEM2-ES, INM-CM4, IPSL-CM5A-MR, and CSIRO-Mk3.6.0. The data is downloaded from the Lawrence Livermore National Laboratory's World Climate Research Programme database powered by the Earth System Grid Federation (Taylor et al., 2012;Cinquini et al., 2014). The five models are selected as they have the required variables in the desired time period and their current or earlier versions have been shown to have better performance relative to other CMIP5 models in simulating the past climate over North America (Radić and Clarke, 2011). The following variables are extracted as daily time series for each GCM: maximum and minimum air temperature, maximum and minimum relative humidity, wind speed, incoming solar radiation at surface, and total precipitation. For the projection period (2080-2099), we retrieved GCM data for two selected emission scenarios, referred to as Representative Concentration Pathways (RCPs; Moss et al., 2010): RCP4.5 and RCP8.5. The different scenarios indicate different amounts of radiative forcing by 2100, where RCP4.5 (RCP8.5) indicates 4.5 W/m 2 (8.5 W/m 2 ). The RCP4.5 emissions scenario is considered a medium stabilization emissions scenario and RCP8.5 is considered a high emissions scenario (Van Vuuren et al., 2011). To obtain the continuous GCM data for 1993-2012 for bias correcting, we use GCM historical runs from 1993 to 2005 and their RCP4.5 and RCP8.5 runs from 2006 to 2012. Evapotranspiration is calculated using GCM data from the Penman-Monteith equations as done with the reanalysis data.
For evapotranspiration, the bias between the present-climate GCM and ERA5 is corrected according to the quantile mismatches within the simulated cumulative distribution function (CDF) using the empirical quantile matching algorithm (Xu, 2015). This method has been shown to outperform other bias correction algorithms when the distribution of a climate variable is unknown or not gamma (Teutschbein and Seibert, 2012;Chen et al., 2013). For precipitation, the bias is corrected using the gamma quantile mismatch between ERA5 and each GCM, with a precipitation threshold of 1 mm, following the method in Xu (2015). Several studies have shown that the gamma quantile matching is one of the better methods for precipitation bias correction (Chen et al., 2013;Lafon et al., 2013).

Water Use and Population Data
To our knowledge, there is no database of daily municipal water use for cities across North America. To collect daily and monthly time series, we directly contacted municipalities. Each time series is at least 3 years in length from the period 1990-2018. We initially collected the data from 55 cities. The datasets from the cities that have strict water consumption laws, large population increases during summer months, or are not continuous or near continuous (missing data with valid observations for previous and subsequent entries) were removed. The gaps in the data are filled with estimates based on linear interpolation, resulting in a set of daily time series of continuous municipal water use for 38 cities and monthly time series of continuous municipal water use for nine cities for a total of 47 cities ( Figure 1B). Though the data has high day-to-day variability which may cause the gap filling to be inaccurate, the gap filling causes little change in the results due to its rareness and the data smoothing applied later. Monthly datasets are converted to daily datasets by using the monthly value for each day within the month.
For each city, we also downloaded its population estimates from the US (United States Census Bureau, 2019) or Canadian (Statistics Canada, 2019a) Census datasets for 1990-2018. Since the datasets provide yearly population estimates, we used a piecewise linear interpolation to derive daily population estimates. County location data is derived from Hauer (2019) and Ramey (2014) for US counties, and Statistics Canada (2019) for Canadian census divisions (Canadian county equivalent). A total of 3,425 counties or county equivalents were used, with 293 coming from Canada and 3,132 from the US.

METHODS
Our overall methodology can be summarized in the following steps: first, we partition the study region into clusters based on similarity in the observed seasonal streamflow and proximity. Second, multiple linear regression (MLR) models, with climate predictors from the ERA5 reanalysis dataset, are calibrated at each river station in order to simulate its seasonal streamflow throughout the 1993-2012 period. Each MLR model is then used to predict the seasonal streamflow regime for 2080-2099 at each station. To do so, the models are forced by stochastically generated downscaled climate scenarios from an ensemble of five GCMs. Third, the same climate data (ERA5) is used to calibrate a set of MLR models to simulate the seasonality in municipal water consumption over the present climate for each of the 47 cities. The locally optimized regression coefficients from the 47 cities are then statistically upscaled to represent the coefficients for each county in the US and Canada. Future changes in seasonal municipal water consumption are then simulated for each county based on the stochastically generated GCM projections. Finally, we aggregate the changes in water supply and demand and assess the risk of water shortages for each region and each season. Seasons are separated into winter (December, January, and February-DJF), spring (March, April, and May-MAM), summer (June, July, and August -JJA), and Fall (September, October, and November-SON).

Clustering
While it is common to partition a large region into smaller domains in order to summarize the results of regional hydrological modelling, there is no consensus on how to define these subregions through clustering (Sawicz et al., 2011). The US Water Resources Council, for example, clustered regions and subregions based on topography, which is the method that has been adopted in assessments of water shortages across the US (e.g., Foti et al., 2012;Mahat et al., 2017). Another common clustering method in hydrological studies is based on climate classification (e.g., Opalinski et al., 2020), or on a combined set of climate and streamflow indices (e.g., Sawicz et al., 2011). As the objective of our paper is to predict seasonality of water shortage and its uncertainty, we chose to partition the study region into clusters on the basis of similarity in streamflow seasonality or regime curves [similar to Brunner et al. (2020)]. For each stream gauge, we derive the average observed streamflow regime over the 20-year current period for each station. The regime curve consists of 366 days in order to account for the measurements during leap years. For non-leap years, the average value from Feb 28 and March 1 is used to fill in for Feb 29. The regime curve of each station is standardized by subtracting it by its mean and dividing it by its standard deviation such that the clustering can focus on regime curve shape similarity rather than on absolute values. We use hierarchical clustering with Ward's algorithm (Ward, 1963), and we determine the optimal number of clusters by trial-and-error in a manner similar to (Knoben et al., 2018). The resulting clusters are then visually inspected on the map in order to further manually partition the clusters whose gauges lie in multiple distinct geographical areas. For example, streamflow regimes in the Pacific and South Atlantic regions were FIGURE 1 | Map of localities with (A) streamflow data (4,290 gauges), and (B) water consumption data (47 cities). In (A) each station dot is colored according to the cluster it belongs to (see Clustering section). In (B) blue and red indicate daily and monthly data, respectively. Some of the city officials asked to keep their cities anonymous, so this map is meant to keep the city names anonymous while showing the spatial extent of the data.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 initially clustered together due to their similar streamflow regimes but were subsequently separated due to their distinct geographical locations. Stations that belong to a streamflow-regime cluster, but are located far outside the cluster's perimeter, are removed as "outliers". Similar to Ivancic and Shaw (2017), we hypothesize that these spatially isolated outliers are present due to significant human impact on the streamflow. Following the steps above, we partitioned the study region into 14 clusters or sub-regions ( Figure 1A). The 4,290 original gauges are reduced to 3,852 due to the removal of outliers.

Streamflow Model
The MLR models are developed to predict the inter-annual variability of streamflow for each gauge using local weather. The MLR models are applied to the averaged streamflow and averaged climate variables (ERA5) using a 30-day moving window. In this way, the model is calibrated across the 20year present period in order to simulate an average value for the same 30-day window (e.g., Jan 1 to Jan 30, Jan 2 to Jan 31, etc . . . ) in each year and for each station. For each region, we validate the model using a validation set of the last 3 years. To obtain the streamflow projections at each station, the model is forced with stochastically generated downscaled climate variables from the ensemble of five GCMs. A well-known limiting assumption of most stationary models is that the model developed over the present climate is assumed valid under future climate conditions (e.g., Ekström et al., 2015). This assumption can be considered acceptable for multiple regression models as long as the variability in present variables (e.g., temperature, precipitation, streamflow) is within the range of their variability in the future period. While the constant variance of the response would be difficult to verify due to the possible non-stationarity and the uncertainty of hydrological conditions over such timescales, we may assume each observation is independent as the observations are separated by a full year such that longer processes including snowmelt and groundwater flow will not affect adjacent years. Below, we provide details on the data smoothing procedure, calibration and validation of the model, and the simulations of future streamflow regimes ( Figure 2).

Calibration and Validation
The MLR models are applied for each of the 366 windows separately, so that the regression is calculated across the 20 (19) points, where each point represents the average value within the current window for the given year. The first model is applied to the first window (1 Jan) across all the years (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), the second model is applied to the second window (2 Jan) across all the years, etc. In this way, we derive FIGURE 2 | Workflow of methods for making predictions for each streamflow station and each window. Boxes indicate data matrices of streamflow (Q), precipitation (P), evapotranspiration (E), or snow days (N). Each row of data represents a single year, and each column indicates a different day.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 366 MLR models in total, each one with its own regression coefficients. The following variables are used as predictors: average daily precipitation (P units in mm), average daily evapotranspiration (ET, units in mm) and total number of days with snowfall (N), where P and ET are derived as the mean over the 30-day moving window, while N is calculated as the sum of days with snowfall >1 mm over the same window. We assume that when evapotranspiration is ≤0, temperature is below or equal to 0°C; consequently, a wet day when evapotranspiration is 0 indicates a snow day. For each window, the modelled streamflow (Q i,j ) (m 3 /s) is derived from the MLR model with optimized regression coefficients (least squares method) as: where a 0,i,j -a 3,i,j are the regression coefficients for station i and window j.
Considering that the MLR models are built on a small sample of 20 points (20 years of data), we use the first 17 years (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)) of data as a training set and 3 years (2010-2012) of data as a test set to test the predictive skill of our methods. We calibrate the model over this reduced sample and predict the streamflow regime curve for the 3 years of left-out data. Similar to Brunner et al. (2020), we then assess the model performance with the Kling-Gupta Efficiency (KGE) between the observed and modelled 3-year streamflow regime curves (Gupta et al., 2009). As the procedure is repeated for each of the streamflow stations, we summarize the model predictive skill by the first quantile, median, and third quantile of the KGEs within each region. Kling has stated that a threshold of 0.5 (0.75) can be used to differentiate poor versus intermediate (good) simulation results using the slightly modified KGE metric, thus we utilize these thresholds to interpret our results (Kling et al., 2012;Thiemig et al., 2013). KGE may be sensitive to outliers and non-normally distributed streamflow (Pool et al., 2018), however, the widespread use of the KGE metric as well as its clear interpretation makes it a strong performance metric.

Postprocessing of Global Climate Model Data
To facilitate the assessment of uncertainties stemming from the GCMs, we input the downscaled precipitation and evapotranspiration time series into a stochastic weather generator so that 2,000 time series of each variable are generated (100 for each of the 20 years). Streamflow for 2080-2099 is then simulated by the MLR models using each of the 2,000 generated time series. For each window, the final simulation is calculated through Monte-Carlo simulations by randomly sampling 1,000 observations from the 100 predicted values from each year to get a large, bootstrapped sample of possible 20-year long time series; then, the values are averaged across the 20 years. In this way, we project 1,000 values of the average seasonal streamflow regime for 2080-2099 for each station and window. Final results are then compressed by representing each day with a mean and standard deviation.
The stochastic weather generator used in this study is developed similarly to the framework illustrated by Wilks and Wilby (1999). The output is generated separately for the two RCP scenarios and for each of the 20 years in 2080-2099. The key steps are summarized as follows: First, for each 30-day window, we have a 5-GCM x 30-day block of data for both evapotranspiration and precipitation. For evapotranspiration, we assume a normal distribution conditioned on the day being wet or dry and derive its mean and standard deviation from the data points for each day type. Second, a transition matrix is calculated from the precipitation data by calculating the four transition probabilities (wet-wet, wet-dry, dry-wet, dry-dry). Then, a two-state Markov chain is created, where p < 1 mm is considered a dry day, and a wet day is assumed otherwise. One hundred 30-day time series of precipitation are generated using the Markov chain to determine wet-dry day occurrences, while the magnitude on wet days is determined by randomly sampling from the precipitation observations greater than 1 mm. The simulation of precipitation magnitude is robust, as no distribution needs to be estimated, so precipitation can be generated non-parametrically. The robust generation of evapotranspiration can be done parametrically because a large set of well-behaved daily data is available for parameter estimation (Semenov et al., 1998;Ababaei, 2014). Finally, evapotranspiration (ET) is derived by randomly generating data from either of the two normal distributions derived earlier dependent on if the Markov chain generated a wet or a dry day. This process generates 100 30-day time series of precipitation and evapotranspiration for each year in 2080-2099 for each 30-day window.

Water Use Model
Daily urban water consumption can be separated into three components that include: long-term (decadal) base water use, calendrical (weekly) water use, and seasonal water use (Wong et al., 2010). The seasonal variability in outdoor water use, which explains a relatively large fraction of variability in municipal water use (DeOreo et al., 2016), can be successfully simulated as a function of weather and climate attributes (Gober et al., 2016;Opalinski et al., 2020). Thus, our working assumption is that the seasonal component of municipal water use, represented by daily time series, needs to be decomposed into its long-term, seasonal, and short-term signal. After the signal decomposition, the model is developed for the seasonal signal only.
While the models are calibrated for selected cities, the objective is to model water use in each county. To do so, the model coefficients are upscaled from city-specific values to values representing each county. Finally, the model is forced by stochastically generated climate scenarios from the five GCMs for both the present and future periods. Below we provide details on the signal decomposition, model calibration, coefficient upscaling, and model projections.

Signal Decomposition
To focus on the seasonal signal in the municipal water use for each city, we remove the calendrical and long-term signals from the city's original time series. Though it may slightly contribute to the violation of observation independence, it is unlikely to distort our results, and removing the calendrical signal is necessary, thus we smooth our time series with a 15-day moving median. This will remove the effects of persistence, holidays, measurement Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 errors, and the day-of-the-week. Long-term variability is mainly driven by population and water efficiency (Wong et al., 2010). To remove the population effects, the smoothed water-use time series are divided by the daily population time series, yielding a daily time series of water efficiency. Each city's water efficiency time series is normalized to be between 0 and 1. A common feature of these time series is their long-term decline that best resembles an exponential decay signal (Supplementary Figure S1). It has been shown that municipal water use per capita in the US has tended to decrease due to government actions that forced manufactures to make waterefficiency improvements for the bathroom and household appliances (Donnelly and Cooley, 2015). These changes were mainly established in the 1990s, so they can be expected to have diminishing returns as time goes on (Donnelly and Cooley, 2015). Empirically, we chose to represent this long-term decline in water use with an exponential decay function that is fitted to the normalized annual water-use time series of each city (Supplementary Figure S1c). Finally, the fitted function, interpolated from annual to daily values, is subtracted from the normalized water efficiency time series (Supplementary Figure S1). The residual daily time series for each city represents the characteristic seasonal municipal water use per capita.

Calibration, Validation, and Upscaling
Each city's remaining time series of daily water consumption per capita represents the response variable in an MLR model, while a set of weather variables from ERA5 data are used as predictors. Prior to training the MLR model, each response variable is standardized (subtracted by its mean and divided by its standard deviation) in order to facilitate the inter-comparison of regression coefficients across cities. After testing different combinations of predictors, we settled with the following list: maximum daily temperature (T, units in°C), daily mean wind speed (U, units in m/s), a binary variable for rain (R, wet/dry day for rain) and a binary variable for snow (W, wet/dry day for snow). Rain and snow days are differentiated using a maximum temperature threshold of 5°C, which is reasonable considering thresholds for mean temperatures are often between −5°C and 5°C (Rajagopal and Harpold, 2016). Standardized per-capita water consumption (C) is modeled with: where b 0,k -b 4,k are the regression coefficients for city k, and R(W) is the binary variable for rainfall (snowfall) which equals 1 if rainfall (snowfall) exceeds 1 mm and 0 otherwise. The model is calibrated for each city, yielding 47 sets of city-specific regression coefficients (Figure 3). Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 Finally, the model coefficients need to be upscaled from their city-specific values to values representing each county. A study with a similar modelling approach to ours found that climatic regression coefficients are correlated in space, with nearby cities yielding more similar regression coefficients than the cities further apart (Opalinski et al., 2020). Thus, we assume that by spatially interpolating the regression coefficients from 47 cities unevenly distributed over the map to the grid of counties, we can obtain a first-order approximation for regression coefficients at each county. We choose to apply smoothed optimal inverse distance weighting with weights set as 1/(1 + d) p , where d is the distance from the county's geographical centre (latitude, longitude) to the city's location, and p is the parameter that needs to be optimized.
We optimized p by maximizing the correlation between interpolated and observed C mod using the following procedure: one out of 47 cities is left out from the sample of known regression coefficients, and its regression coefficients are instead derived from the interpolation method with each possible p. This procedure is iterated until all 47 cities are removed once from the original sample. The final value of p is the one that yields the largest mean correlation between interpolated and observed C across all the cities. The optimization yielded p 2.1 (Supplementary Figure S2).

Post-processing of Climate Data
Time series of normalized water consumption per capita (C) for each county for both the present (1993-2012) and future (2080)(2081)(2082)(2083)(2084)(2085)(2086)(2087)(2088)(2089)(2090)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098)(2099) periods are simulated by the MLR model forced by stochastically generated climate data. Similar to the supply model, data from the five GCMs and two RCP scenarios are used as inputs to the stochastic weather generator. Present climate data is patched from the GCM historical runs from 1993 to 2005 and their RCP4.5 runs from 2006 to 2012. Because changes in the water use per county are assessed as a difference between modelled future and modelled present water use, from the same ensemble of GCMs, no bias-correcting is applied.
The predictors in the model (Equation 2) are stochastically generated in the following way: precipitation wet/dry days are derived from a two-state Markov chain where the transition probabilities are estimated from a 31-day moving window of GCM precipitation time series. For temperature, using the same 31-day moving window, we assume a normal distribution and derive its means and standard deviations conditioned on the presence of precipitation. For wind speed, a Weibull distribution is chosen following previous studies (e.g., Seguro and Lambert, 2000;Alizadeh et al., 2019), and we also assume that the wind speed is conditioned on the day being wet or dry. As a result of this assumption, two sets of two-parameter Weibull distributions are estimated using maximum likelihood estimation. The differentiation between snowfall and rainfall is determined by a maximum temperature threshold of 5°C as was done in the calibration phase.
One hundred time series of each weather variable is generated for each year (100 × 20 sets of 365-day time series). Each of the 2,000 generated sets of weather time series are then used to predict a 365-day time series of C. To obtain the uncertainty estimate for C for each day, we randomly select 1,000 outcomes from each of the 100 time series of each year. Like with the watersupply model, we take the average C across the 20 years, so a total of 1,000 20-year average time series of normalized water consumption per capita for each county are simulated. The final results are presented as the mean ± a standard deviation of C for each day in the present period, and the mean ± a standard deviation of C for each day in the future period for both RCP4.5 and RCP8.5.

Risk Scores
Defining vulnerability or risk of water shortages is difficult to do so precisely while also being implementable (Foti et al., 2012). Due to this limitation, each previous assessment of water shortage risk has developed their own risk assessment score (Hurd et al., 1999;Foti et al., 2012;Roy et al., 2012;Dickson and Dzombak, 2019). For example, Foti et al. (2012) defined vulnerability as the probability that demand will be less than supply in a region, while Dickson and Dzombak (2019) defines an overall risk index from five components of risk including 1) annual proportion of use of local water supply, 2) summer proportion of use of local water supply, 3) projected increase in demand, 4) summer proportion of use of local water demand, and 5) proportion of groundwater withdrawal to total withdrawal.
Building on previously proposed risk scores, we introduce a metric as a proxy for the risk imposed by climate change on water resources for each subregion in the study domain. The metric, defined as the risk score (RS), combines the projected changes in water supply for each station with the projected changes in water demand for each county, and congregates them by sub-region in the following way: where D is the probability of a water demand increase from 1993-2012 to 2080-2099, CD is the normalized demand, S is the probability of a supply decrease from 1993-2012 to 2080-2099, and SC is one minus the normalized supply. Each river station is already assigned a sub-region (Figure 1), while each county is assigned to the sub-region to which the closest river station belongs. The probability D is calculated by forming a new normal random variable, X, for each RCP scenario with the mean determined by the projected demand in 2080-2099 subtracted by the estimated demand in 1993-2012 and variance set as the sum of the variances of projected demand and estimated current demand. S is calculated in the same way as D, but using the projected streamflow in 2080-2099 and the streamflow regime curve from 1993-2012 to form a random variable Y. In summary, for day i, region j, county k and streamflow station l: D i,j mean k P X i,j,k ≥ 0 (4) S i,j mean l P Y i,j,l ≤ 0 D and S become a 365-day x 14-region matrix with values between 0 and 1. To compute CD and SC, the demand and supply regime curves for 1993-2012 are scaled to be between 0 and 1 by Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 subtracting by their minimums and dividing by the difference of their maximums and minimums. For each subregion, final CD is calculated as the mean CD across all stations belonging to the region, while final SC is calculated as the mean SC across all counties belonging to the region. Thus, each final CD and SC is a 365-day x 14-region matrix with values between 0 and 1. RS is thus a 365-day x 14-region matrix with values between 0 and 4, where 0 implies no risk and 4 implies maximum risk of water shortages.

Regional Clusters
The clustering of streamflow data across the region yielded 14 clusters or subregions (Figure 4). Below we describe the main characteristics of the average streamflow regime and some general climatic attributes for each cluster: 1) The first cluster, which we name "the West Coast", is characterized by high winter streamflow and minimum flows in late August. The cluster includes Vancouver Island, western half of Washington, Oregon, and California. The climate in these areas is characterized by high frequency of precipitation events in winter that are predominately rainfall rather than snowfall. Thus, relative to other clusters, the streamflow is mainly driven by rainfall. 2) The "Southwest" cluster is comprised of the desert areas of California and Arizona. Although very little precipitation is available in these areas, precipitation may occur periodically in winter months. These rainfall events are usually of high intensity, potentially causing infiltrationexcess overland flow. Thus, these catchments are rainfall dominated and have relatively quick responses to precipitation input. 3) The "Deep South" cluster includes Louisiana, Mississippi, Alabama, Georgia, Tennessee, and South Carolina. These states are known to have warm and humid summers with some precipitation. Winter and spring have the highest precipitation. The landscape is covered by forests which implies overland flow is less likely compared to the Southwest cluster. With such active storage, the streams' Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 response to precipitation may be slower despite being a raindominated region with a peak in the late winter months and low flows in late-summer. 4) The "Mid-Atlantic" cluster includes North Carolina, Kentucky, Virginia, West Virginia, Ohio, Pennsylvania, Maryland, Delaware, New Jersey, southern New York, Connecticut, and Southern Ontario. Mostly covered in forest, this region is known for its hot summers and very cold winters, where winter temperatures are often at or below freezing. Precipitation is often high, with precipitation being fairly evenly distributed throughout the year. This region may be hydrologically dominated by rain-on-snow events in winter, with all snow melting by early spring. The streamflow regime curve in this region has a peak in early spring and a low in mid-summer. 5) The "South-Central US" consists of Northeast Texas, Southeastern Oklahoma, Arkansas, and northern Louisiana. These areas are covered by grass, shrubs, with some forested areas further east. Precipitation is usually highest in the spring in these areas, with little to no snowfall, leading to highest streamflow in the spring with a seasonal drought in mid-summer. 6) The "Florida" cluster is largely covered by wetlands in the southern half and forest in the north. Since Florida is subtropical and tropical in some areas, precipitation occurs throughout the year. Summers are slightly wetter and more humid than are winters. The regime curve in this region has little variance, with highs in the late summer and lows in early summer. 7) The "Hawaii" cluster is in the tropics with precipitation and warm, humid weather throughout the year. The streamflow regime curve in this region has extremely low variability, with no clear high or low. 8) The "Rockies" cluster consists of Northern Wyoming, Idaho, Montana, Eastern Washington, and parts of British Columbia and the Yukon. This region is extremely cold in winters, often with extensive snow cover. The streamflow stations are in snow-dominated catchments, so the hydrographs are characterized by low streamflow in winter when no snowmelt is present, followed by peak streamflow in late spring to early summer as snowmelt intensifies. 9) The "Alaska" cluster, similarly to the Rockies, is a snowdominated region; however, in Alaska, the colder climate and extensive snow cover, as well as glacierized terrain, allows for the snowmelt and/or glacier runoff to last throughout the summer; thus, giving hydrographs with a wide peak throughout the summer. 10) The "Midwest" cluster contains Illinois, Wisconsin, Iowa, Minnesota, North and South Dakota, as well as Southern Alberta and Saskatchewan. This region's land use mainly consists of agricultural areas and grasslands. This region is fairly dry; nevertheless, snow cover persists throughout much of the winter, with melt occurring in the spring leading to the peak in streamflow. 11) The "Northeast" cluster consists of New York, Vermont, New Hampshire, Maine, Southern Quebec, New Brunswick, Nova Scotia, Prince Edward Island, and Southern Newfoundland. This region is covered by forest and has very cold winters; hence, snow builds in the winter and then snowmelt generates peak streamflow in April. 12) The "West" cluster includes parts of British Columbia, Alberta, Idaho, Montana, Utah, Nevada, California, and Arizona. The hydrographs in this region peak in spring as a result of snowmelt. 13) The "Central US" cluster consists of Southern Minnesota, Iowa, Nebraska, Missouri, Kansas, and Oklahoma. The hydrographs in this region peak in spring. 14) The "Eastern Canada" cluster includes northern Newfoundland, Quebec, Ontario, and Manitoba. The hydrographs in this region peak in spring. Winters are cold and snow cover is prevalent, so runoff is driven by spring snowmelt.
According to these characteristics we can group the 14 subregions into three general regime categories: rainfall dominated streamflow regime (West Coast, Southwest, Deep South, South-Central US, Florida, and Hawaii), rain-on-snow dominated regime (Mid-Atlantic, Midwest, and Central US), and snow dominated regions (Rockies, Alaska, Northeast, West, and Eastern Canada).

Model Performance
The KGE is estimated for all streamflow stations within each region for the 3-year test set ( Table 2)

Projections of Streamflow
An example of modelled future (2080-2099) versus present (1993-2012) streamflow regime curves are shown in Figure 5 for a streamflow station in California. The streamflow at this station is expected to substantially increase in late winter and early spring, while the streamflow will remain almost unchanged throughout the rest of the year. The uncertainty interval for both RCP scenarios is large, especially for the time window with the largest projected increase in streamflow. We estimate the watersupply components of the risk score for this station (Equation 3): on a seasonal scale, the largest CS is in the fall (median CS 1.00) and the smallest is in the winter (CS 0.26), while the likely decrease in streamflow is projected for the summer (S 1.00 for RCP4.5; S 1.00 for RCP8.5), and the least likely decrease is projected for the winter (S 0.27 for RCP4.5; S 0.27 for RCP8.5).
Next, we summarize the results of regional streamflow projections for each season, with more focus on subregions where the model has higher predictive capabilities (larger KGE) relative to other subregions: Results for winter streamflow (DJF; Figure 6) reveal that the likelihood of a downward change in streamflow is linked to latitude. Specifically, for low-latitude regions (Hawaii, Southwest, Florida, Central US), streamflow is projected to decrease for one or both RCP scenarios. For these regions, the average S is 0.56 (0.68) for RCP4.5 (RCP8.5) scenario, while the average CS is 0.70. Out of these regions, only Central US and Florida have KGE <0.5 implying a low confidence in the projections. The probability of a decrease in Hawaii is certain with the interquartile range (IQR) of 0 for both RCP4.5 and RCP8.5 scenarios. The IQRs for Southwest, Florida, and Central US are much larger, ranging from 0.35 for Southwest's RCP8.5 scenario to 0.95 for Florida's RCP4.5 scenario. Streamflow is projected to increase in Alaska, the Rockies, the Midwest, the West Coast, and the Northeast. According the KGE values, however, we have higher confidence in the projections for the West Coast (KGE 0.82), Alaska (KGE 0.89) and less certain for Rockies (KGE 0.73), Midwest (KGE 0.59), and Northeast (KGE 0.70).
In simulating the streamflow in spring (MAM), under the RCP4.5 scenario, Hawaii and the West Coast will experience decreases in streamflow (their average S is 0.76 and there IQRs are 0 and 0.48), whereas the Deep South, East Coast, Southern Central US, Central US, Northeast, and the Midwest will experience increases (their average S is 0.13 and their IQRs are 0.19, 0.11, 0.25, 0.17, 0.38, and 0.10 respectively). Like the RCP4.5 scenario, RCP8.5 will lead to decreases in streamflow in Hawaii, West Coast, Florida, and the Southwest (their average S is 0.75 and their IQRs are 0, 0.59, 0.12, and 0.38 respectively). These results are particularly of note for Hawaii, the Southwest, and Florida as these regions experience low present water supply in the spring relative to other seasons (their average CS is 0.64).
Projections of summer (JJA) streamflow reveal that water resources will experience strain across most of North America ( Figure 6). According to the RCP4.5 scenario, the West Coast, the Southwest, South Central US, Florida, the Deep South, Central US, and the Midwest will all likely experience decreased water supply (their average S is 0.80). Following the RCP8.5 scenario, in addition to the regions affected by the RCP4.5 scenario, Hawaii, the East Coast, the Northeast, and Eastern Canada will see decreases in summer streamflow (their average S is 0.81). The IQRs of decrease probabilities are small for the West Coast (0 for RCP4.5, 0 for RCP8.5), Southwest (0.18, 0), Florida (0.08, 0), and Central US (0.18, 0).
For the projections of autumn (SON) streamflow, regardless of the RCP scenario followed, eastern regions are most likely going to be affected by decreases in their streamflow, including the Deep South, the East Coast, Southern Central US, Florida, Midwest, Northeast, and Central US (their average S is 0.73 and 0.69 for RCP4.5 and RCP8.5, respectively). The IQRs of decrease probabilities are small for the Southern Central US (0.05 for RCP4.5, 0 for RCP8.5) and the Deep South (0.20, 0.10). Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631

Simulations of Water Demand
Model Performance Supplementary Figure S2 shows the mean correlation (0.69) between predicted and observed daily water efficiency with optimal inverse distance weighting calculated using leave-oneout cross validation. All locations have at least 1,000 observations of daily water use efficiency, so this correlation is significant (average p-value < 0.01), and the predictions are representative of the observations.

Projections of Municipal Water Use
An example of modelled future versus present water consumption per capita is shown in Supplementary Figure  S3 for a county in Alabama. As illustrated in the figure, there is a pronounced seasonal cycle in the water consumption, with the highest consumption during the summer months. The changes in water demand are projected to occur in all seasons, with the largest increases in summer, and the lowest in winter. These results are reflected in the water-use components of the risk score (Equation 3): the largest CD is in summer (CD 0.95) and the smallest CD is winter (CD 0.08).
The likelihood that this county will experience an increase in demand is also the largest in summer (D 0.97 for RCP4.5; D 1.00 for RCP8.5), and the smallest in winter (D 0.86 for RCP4.5; D 0.98 for RCP8.5). Projections for each subregion and each season reveal that the changes in water consumption are relatively uniform across seasons, with the main difference being between winter and summer water consumption (Figure 7). For all regions, the current water demand is the largest in summer and lowest in winter (their average CD is winter is 0.07, and their average CD in summer is 0.92). In all regions, RCP4.5 will substantially increase the annual water demand (average D 0.93 for RCP4.5), while more increase is expected for RCP8.5 (average D 0.99 for RCP8.5). Most of this increase is expected in the summer season (average summer D 0.95 for RCP4.5; D 1.00 for RCP8.5). The top three subregions with the highest projected increase in demand in summer given the RCP4.5 and RCP8.5 scenarios are Hawaii and Florida, and their average summer D is 0.99 and 1.00, respectively. The top three subregions with the highest projected increase in demand for RCP4.5 (RCP8.5) in winter are: Hawaii, Southwest, and West Coast (Hawaii, Southwest, and Eastern Canada), and their average winter D 0.96 (1.00). During all seasons and in all regions, increases in municipal water demand are expected, regardless of the RCP scenario followed, however the increases are more certain under RCP8.5 compared to RCP4.5. The interquartile range of all regions and seasons given the RCP4.5 scenario ranges from 0 to 0.1, whereas the IQRs given RCP8.5 remain at 0.

Risk Assessment
Risk scores (RS) for each region and RCP scenario are quantified for each season. Because of consistently higher demand in summer (Figure 7), total risk is consistently highest in the summer, so only summer risk scores are shown (Figure 8).
According to the RCP4.5 scenario, the region with the highest risk of freshwater shortage in the summer (highest median RS) is the Southwest (region 2), followed by the West Coast (region 1) and the Deep South (region 3). The West Coast and Southwest both have low uncertainty (IQRs of 0.30 and 0.10), but the Deep South has high uncertainty (IQR 0.55) stemming from high uncertainty of streamflow projections in this region. According to the RCP8.5 scenario, the regions with the highest risk are: Southern Central United States (region 5), the Deep South (region 3), and Northeast (region 11). The IQRs of these risk scores are low for Southern Central US (0.40) and the Deep South (0.37).

Summary of Results
GCMs generally agree that the entirety of North America will experience temperature increases throughout this century, while precipitation will change over time as a function of latitude, where southern areas will see decreases in precipitation, and northern areas will experience increases (Roy et al., 2012;Forster et al., 2020). This will likely lead to decreases in the fraction of precipitation falling as snow, increases in evapotranspiration, and changes in available water. But how this change in climate will alter the supply and demand of freshwater is poorly understood. Here we aimed at building a better understanding of these uncertainties while obtaining a general overview of where and when the risk of supply deficits is greatest. We focused our paper on two objectives: First, we aimed to determine the viability of using simple regression models to capture the variability in streamflow and water consumption using globally available climate reanalysis data while assessing the risk of water shortage in multiple regions across North America. The viability of simple regression models for the use of modeling streamflow varied by region, as measured by the validation phase KGE. Florida and Central US had the lowest median KGEs. The other regions had highly acceptable KGEs, especially when the regime curves do not change much from year to year, as is the case in Alaska. Further, the correlation between simulated and observed water use was high. This result means that urban water use in a city can robustly be estimated from population data and climate data in conjunction with water use data from surrounding cities. We found that demand is consistently highest in the summer, thus making water shortage risks consistently highest in the summer. In summer, the highest risk scores are in the Southwest (Southern Central US), West Coast (Deep South), and Deep South (Northeast) given the RCP4.5 (RCP8.5) scenario.
As a second objective, we sought to quantify the uncertainties in predicting streamflow and water consumption and therefore water shortage stemming from GCM climate projections. We found that almost all of the risk score uncertainty is derived from the uncertainty of streamflow projections from GCM data, while there is little-to-no uncertainty in demand projections from GCM data. This occurs because the supply is driven by precipitation while demand is driven by maximum temperature, and GCM projections of temperature have a narrower spread of projections from the weather generator compared to precipitation.

Comparison With Previous Works
Similar to previous studies that used regression models for predicting streamflow (e.g., Barbarossa et al., 2017;Mendoza et al., 2017), we found the MLR models to be reliable in most of our 14 subregions. However, in some regions (e.g., Florida and Central US), the model performance was relatively poor. The poor performance in these regions is not limited to simple regression models. Brunner et al. (2020), using physically based models, also found that Central US was the most difficult region to reliably simulate water supply.
For the US, our projections of streamflow generally agree with the previous findings that the annual streamflow will substantially decrease in South and Western US (Hurd et al., 1999;Foti et al., 2012;Chien and Knouft, 2013;Mahat et al., 2017). Throughout the region, summer streamflow was shown to be more susceptible to decreases (Maurer and Duffy, 2005;Chien and Knouft, 2013) than streamflow in other seasons. In particular, under the RCP4.5 scenario, we projected that Florida may experience significant streamflow decreases in the summer, corroborating previous findings in (Mahat et al., 2017). Furthermore, under the RCP8.5 scenario, our results also agree with Mahat et al. (2017) that the Southwest, Southern Central US, and Florida will experience substantial drops in streamflow. These projections are relatively consistent across the climate projections stochastically generated from the five GCMs.
Our projections of water use agree with those is Foti et al. (2012), showing that climate-driven water use will increase across the United States, as well as Canada, and uncertainties around these projections will decrease as temperatures rise further. Since maximum temperature is the dominant driver of water consumption (Opalinski et al., 2020), and all CMIP5 models Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730631 agree that temperatures are projected to rise across the entire United States (Mahat et al., 2017), this is in-line with expectations. The MLR modeling of municipal water demand showed that most of our city-specific regression coefficients are highly spatially dependent as previously found in (Opalinski et al., 2020). The regression coefficient that shows the contribution of the number of snowy days to water demand is negatively correlated with latitude, implying that the impacts of snowy days on water demand is greater in lower latitude areas. For example, at low latitudes, the snowfall-day coefficient is between 0.6 and 1.4, while in Canada, the coefficient drops to 0.2. The regression coefficient that shows the contribution of the maximum temperature is also spatially varying, with the highest values in western and southern areas (0.08-0.19), and the lowestin northern and eastern areas (0.03-0.08). In the West Coast, southwest US, and Florida regions, the coefficient multiplying the wind speed is normally positive (0-0.15), while in other places, the coefficient of wind is normally negative (−0.15-0), implying that warm winds could raise water demand, while colder winds could decrease demand. However, the coefficient in Texas is negative, so this hypothesis should be tested further. Rain-day coefficients seem to not be spatially dependent. We compared our spatially distributed coefficient values with those derived in Opalinski et al. (2020) and found general agreement for the coefficient multiplying the maximum temperature. In Opalinski et al. (2020), however, the link with the latitude is not as clear as in our case which is likely because we extended the assessment to Canadian cities while their study only provides the results for the contiguous US. Regarding precipitation, Opalinski et al. (2020) found that the increase in precipitation in cold regions was associated with an increase in water consumption. One interpretation for this finding is that the increased snowfall is linked to increased water use, which is also consistent with our results. To our knowledge, there have been relatively few up-to-date assessments of water resource vulnerabilities and shortage risks due to climate change across the entire US or North America. Dickson and Dzombak (2019) as well as Foti et al. (2012) used physically based models to project streamflow, while Roy et al. (2012) calculated available precipitation (precipitation minus evaporation) instead of streamflow. Previously, urban waterdemand projections have been made based on changes in population and electricity production (Roy et al., 2012;Dickson and Dzombak, 2019), or climate regression techniques combined with the effects of income, population, and efficiency (Foti et al., 2012). Our study is the first that focused on both the US and Canada, and the first to combine the use of regression models for both supply and demand projections.
Our results generally agree with those in the previous studies for the regions where the previous projections of risks have been assessed. For example, the probability of water shortage was calculated across the U.S. by Foti et al. (2012) and probabilities were highest in the same areas as our study indicated, namely the southwestern United States and the coast of California, though they showed a high risk for the central US as well. Roy et al. (2012) found that the southwest and central US as well as Florida will have the highest risk index in 2050. The latest assessment of risk showed that central California, and parts of Nevada have the highest risk of water shortage (Dickson and Dzombak, 2019). Regardless of the methods, data, or scenario chosen, all of these studies agree that California and the US Southwest have the most troubling futures, relative to other parts of North America, when it comes to the effects of climate change on water resources.

Limitations
Though we used an expansive collection of streamflow data, continuous, long-term hydrologic data is often limited to large rivers and our conclusions may not accurately extrapolate to smaller basins (Kovach et al., 2019). In a similar fashion, the water use data is mostly limited to medium-sized cities and the captured relationships between water use and weather may not extrapolate to extremely large cities or small towns. Although some hydrological extremes are included in the 20-year period of analysis used in this paper (e.g., the 2000-2004 drought in western North America), the length of time used may be seen as too modest as the full extent of climate extremes in all areas may not be present and therefore learnable in our data (Schwalm et al., 2012). Rather than developing the best possible prediction model, we aimed at exploring how well a relatively simple model (MLR) can explain interannual and seasonal variability in water consumption and streamflow. Likewise, our water shortage vulnerability framework is relatively narrow in scope. For a more complete risk assessment, one would have to consider biological and ecological sensitivity as well as adaptive capacity (Kovach et al., 2019). Though we utilized stochastic weather generators to increase the number of cases in our ensemble of future scenarios, we only included five GCMs as a baseline, so the results could be sensitive to our choice in GCMs.
Due to the number of models included for projecting streamflow, the assumptions of linear behavior and constant variance over the long period of study could not be verified. The least success in terms of MLR performance among the streamflow models is found for Central US and Florida. One of the key reasons for the poor performance is the fact that runoff could be controlled by non-linear threshold-like belowground hydrologic processes (Jones et al., 2019). In some catchments, precipitation can fall and either quickly flood the surface when the soil is saturated or recharge the groundwater system. Thus, water from similar sized precipitation events may be transferred toward the stream immediately or a few months to several decades later. The "memory" in the hydrological system, which acts on multiple time scales, is not possible to be captured by the linear regression model used in our paper. In deforested, arid, urbanized, or agricultural areas with little snowfall and little interactions between precipitation and belowground attributes, streamflow response to precipitation is relatively fast. In these regions, as well as in regions with similar year-to-year regimes, our ensemble of regression models is deemed acceptable to forecast future water shortage and its uncertainty.
With regards to our projections of water use, unlike previous studies, we did not consider water use from sectors such as agriculture or thermoelectric power which admittedly take up a large portion of water (Foti et al., 2012;Roy et al., 2012;Dickson and Dzombak, 2019). Climate successfully explained a large portion of the variability in our urban water use data, and climate change will likely explain much of the future changes in use, however, non-climatic changes can also be substantial factors. For instance, age of the population, household size, and other sociodemographic attributes can explain some water use behavior (March and Saurí, 2010). While the effects of population, income, and technological efficiency growth have been shown to negate overall, individual regions or cities may experience outsized effects of one or another (Rosegrant and Cai, 2002). City or building specific attributes related to water use such as building type, tree fraction, building age, impervious surface percentage, and population density can also change over time and should be seriously considered (Stoker and Rothfeder, 2014;Chang et al., 2017). General attitudes towards conservation and lawn upkeep are yet other examples of features that may change in tandem with water use (Hong and Chang, 2014). Further, we assumed that spatially interpolating our estimated regression coefficients from cities with data to the grid of counties gives a good enough first-order approximation of the effects of climate on water use in each county. Though the validation provides some confidence that interpolating the regression coefficients yielded optimal results for the purpose of our study, some subregions had little-to-no data (e.g., Hawaii and the Deep South), and therefore uncertainty in their results remains unaccounted for.

CONCLUSION
Here we presented a new method for quantifying future risks of municipal water shortages across North America in response to climate change. We applied a set of multiple regression models to project the changes in streamflow regimes (water-supply) and urban water use regimes (water-demand) in response to climate scenarios stochastically generated from an ensemble of five GCMs, for the RCP4.5 and RCP8.5 emission scenarios. The models were calibrated over the present period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and then used for projections over the 2080-2099 period. The results were analyzed for each of 14 identified clusters or subregions across the original domain (the US and Canada).
Results show that risk scores are considerably higher in the summer compared to the rest of the year, since urban water demand is high regardless of the location chosen. The resulting risk scores for water shortages by the end of the 21st century under the RCP4.5 emissions scenario, indicate a high risk for the West Coast, the Deep South, and the Southwestern US. Under the RCP8.5 scenario, the regions with the highest risk are Southern Central US, the Deep South, and Northeast, while the regions with the lowest risk are Alaska, Rockies, and Hawaii. The high risk scores are due to high and increasing demand concurrently happening with low and decreasing supply. Almost all uncertainty in the risk scores is rooted in the uncertainty of streamflow projections from GCM data, while there is little-to-no uncertainty in water-demand projections from GCM data. Overall, the predictive power of our streamflow model is shown to be sufficient in most subregions, while for Florida and Central US, the model yielded low predictive power. Our study reveals that simple statistical models can produce projections of water shortages that agree with those derived from more complex methods used in previous studies. Nevertheless, more work is needed to identify and implement mitigation strategies to prevent water resource shortages. In particular, the regions that we identified as those with high risk scores, should be examined in more detail so that their water resources can be analyzed at the community level. Future research should also aim to employ more complex methods, especially in those regions shown with our simple method to have high risk and relatively large uncertainty in predicted water shortages. To capture the long-term 'memory' and non-linear threshold-like behavior of hydrologic systems, due to snowmelt characteristics or differing geological water partitioning systems, more complex and computationally expensive methods are needed, such as Long-Short-Term Memory models (e.g., Shen, 2018), as well as diverse climatic and physical data.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: CMIP5 data is available via https://esgf-node. llnl.gov/search/cmip5/. ERA5 data is available via https://cds. climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-singlelevels?tab overview. Streamflow data for US and Canada can be obtained from USGS (https://waterdata.usgs.gov/nwis/sw) and the HYDAT database (https://www.canada.ca/en/environmentclimate-change/services/water-overview/quantity/monitoring/ survey/data-products-services/national-archive-hydat.html). The water use data for each city has been provided by individual municipalities and is not currently publicly available, however, some of the data can be shared if the corresponding author is contacted directly.

AUTHOR CONTRIBUTIONS
JJ carried out data collection and data analysis. JJ and VR drafted the manuscript. JJ, VR, and AA designed the methods and edited the manuscript.