## Abstract

A gridded hydrologic model was developed to assess the impact of projected climate change on future Delaware River Basin (DRB) hydrology. The DRB serves as a water supply resource to over 15 million people. Model evaluation statistics for both water year and monthly runoff projections indicate that the model is able to capture well the hydrologic conditions of the DRB. Basinwide, annual temperature is projected to increase from 2.0 to 5.5 °C by 2080–2099. Correspondingly, potential and actual evapotranspiration, precipitation, rainfall, and runoff are all projected to increase, while snowfall, snow water storage, snowmelt, and subsurface moisture are all projected to decrease. By 2080–2099, basinwide summer subsurface moisture is projected to decrease 7–18% due to increased evapotranspiration, while winter runoff is projected to increase 15–43% due to increased precipitation and snowmelt and a conversion of snowfall to rainfall. Significant spatial variability in future changes to hydrologic parameters exists across the DRB. Changes in the timing and amount of future runoff and other hydrologic conditions need to be considered for future water resource management.

## INTRODUCTION

The Delaware River Basin (DRB) is a 35,066 km2 watershed that contains parts of the states of New York, Pennsylvania, New Jersey, and Delaware along the mid-Atlantic coast of the USA (Figure 1). The DRB is a drinking, agricultural, and industrial water supply to over 15 million people (approximately 5% of the US population) (DRBC 2019a). Nearly half of these people live outside of the basin, including the residents of New York City who receive almost half of their water supply from three reservoirs on tributaries of the Delaware River (DRBC 2019a). The main stem of the Delaware River, at 531 km, is the longest undammed river in the United States east of the Mississippi River. Three sections of the main stem and many sections of the 216 tributaries that feed the Delaware River are part of the National Wild and Scenic Rivers System. Economic activity throughout the watershed generates approximately $22.5 billion per year, while activity in the ports along the Delaware Bay at the outlet to the Atlantic Ocean generates$19 billion per year (DRBC 2019a).

Figure 1

Maps showing (a) political boundaries, (b) elevation, (c) available water storage, (d) subwatersheds for validation, and (e) 2000 LULC for the modeled domain. For (d), gages 8 and 10 are highlighted as they are referred to in following figures. For (e), wat, water; dev, developed; for, forest; agr, agriculture.

Figure 1

Maps showing (a) political boundaries, (b) elevation, (c) available water storage, (d) subwatersheds for validation, and (e) 2000 LULC for the modeled domain. For (d), gages 8 and 10 are highlighted as they are referred to in following figures. For (e), wat, water; dev, developed; for, forest; agr, agriculture.

Due to its status as an economic engine and a water supply entity for the region, the goal of this study is to examine the impacts of projected future climate change on the hydroclimatology of the DRB. Several government and nonprofit organizations have summarized the historic and projected future changes in climate for the DRB (Beecher et al. 2013; Andrus et al. 2017; Haaf et al. 2017). Over the last century, the annual mean temperature for the DRB has increased 1.0–1.2 °C. The trend is greater in the recent past when compared to the entire time series. Warming is significant across all seasons. By the end of the 21st century, temperatures are projected to warm 1–5 °C with the greatest warming in winter and secondarily in summer.

Historic trends in annual precipitation totals, while slightly increasing, are not statistically significant for the DRB. These trends are also larger when comparing the recent past to the entire time series. Seasonally, only increases in fall precipitation are statistically significant. In the future, precipitation is projected to increase 10–14% in winter, with lesser increases in other seasons.

Numerous studies have examined the impact of historic and projected temperature and precipitation changes on the hydrology of the mid-Atlantic region (Najjar 1999; Najjar et al. 2000, 2009, 2010; Neff et al. 2000; Smith & Hawkins 2011; Hawkins & Austin 2012; Bastola 2013; Hawkins 2015) and the DRB (McCabe & Ayers 1989; McCabe & Wolock 1992; Neff et al. 2000; Frei et al. 2002; Burns et al. 2007; Najjar et al. 2009; Matonse et al. 2011, 2013). Earlier studies (McCabe & Ayers 1989; McCabe & Wolock 1992; Neff et al. 2000; Frei et al. 2002; Najjar et al. 2009) suggested a potential decrease in the DRB and/or Catskill Mountains (location of New York City's primary water supply) runoff due to increasing temperatures and great uncertainty in precipitation projections. From a water supply standpoint, relatively large increases in precipitation would be needed to offset increased evaporation due to warming. More recently, Burns et al. (2007) looked at Catskill Mountain hydroclimate trends from 1952 to 2005 which showed increases in temperature, precipitation, potential evapotranspiration (PE), and runoff, although PE values are likely overestimated (Milly & Dunne 2017). Using improved global climate model data, Matonse et al. (2011, 2013) showed that in the future, increased winter temperature, increased winter precipitation, and earlier snowmelt all combine to increase winter runoff and decrease spring runoff.

Most recently and most relevant to the current project, Williamson et al. (2016) used the WATER model (Williamson et al. 2015) to assess future hydrologic changes in three tributaries of the DRB and the impact of different formulations of PE in those assessments. The simulations showed increased annual runoff, resulting mostly from winter runoff increases due to higher winter precipitation and a shift away from storing water as snowpack in the winter. Importantly, the choice of PE formulations had a dramatic impact on runoff results. The energy-based (Priestley & Taylor 1972) formulation produced annual increases in runoff, while the temperature-based (Hamon 1961) formulation produced smaller increases or even decreases in runoff. Thus, parameterization of PE is a critical component of assessing future hydrologic change.

While the previous studies have examined the historic and future hydrology of the DRB, they differ from the present study in a few significant ways that are presented below and will be elaborated upon in the Methods section. Consequently, the resulting data have great potential utility for DRB stakeholders.

• 1.

The current hydrologic model runs on a distributed grid for the entire DRB where previous studies have examined the entire DRB as a whole or selected sub-basins. Therefore, geographic variability for the DRB may be ascertained.

• 2.

The current hydrologic model is driven by 231 future climate projections based on 37 different climate models, 4 different representative concentration pathways (RCPs), and several runs for different model/RCP combinations where previous studies have used only a few climate models and one or two RCPs.

• 3.

The current hydrologic model produces gridded time series of all hydrologic data for all 231 climate projections where previous studies have only calculated hydrologic values for selected future time slices and for coarser resolution geography.

Given the economic and water supply reliance on the DRB, it is critical to understand how future climate change may impact the resource. At both the national (EPA 2011, 2016) and local DRB level (Matonse et al. 2013; DRBC 2017), agencies are planning for how future climate change will impact freshwater supplies. Considerations include but are not limited to reservoir supply (Matonse et al. 2013), drought (Kang & Sridhar 2017), and future climate extremes (Ning et al. 2015). The human impacts of these future changes and water resource management considerations can only be properly understood with accurate assessments of these resources.

## METHODS

### Study area and model physical parameters

Figure 1 shows the hydrologic model domain, which was selected as the area that encompassed all the counties that intersected with the DRB (Figure 1(a)). The model was gridded based on the resolution of the input climate data (discussed below). Grid cell resolution was 1/8° latitude by 1/8° longitude (∼12 km). The entire model domain was 36 rows × 24 columns and contained 864 grid cells. The DRB itself contained 225 grid cells.

Elevation data at a 30 m resolution were obtained from the National Hydrography Dataset (USGS 2015) and aggregated to the 12 km resolution associated with the hydrologic model (Figure 1(b)) by averaging the 30 m pixels. Elevations range from the sea level along the coastal plain and the Delaware Bay in the south to 930 m in the Catskill Mountains in the northern part of the DRB. Subsurface moisture storage for the hydrologic model was calculated from the available water storage (AWS) field in the Gridded Soil Survey Geographic database (gSSURGO; USDA 2017). gSSURGO data originally came at a resolution of 10 m and were aggregated to 12 km (Figure 1(c)) by averaging the 10 m pixels. Subwatersheds to evaluate modeled streamflow were selected to correspond to major tributaries of the Delaware River, to represent a range of latitudes, elevations, sizes, and physiographic regions, and to align with subwatersheds used in other components of this broader project (Woltemade 2016; Woltemade et al. in review) (Figure 1(d) and Table 1).

Table 1

Description of stream gages in Figure 1(d) and the entire DRB

No.USGS IDNameDrainage area (km2)Drainage area ratioa
01413500 East Branch Delaware River at Margaretville, NY 422.2 0.737
01415000 Tremper Kill near Andes, NY 86.0 0.601
01445500 Pequest River at Pequest, NJ 274.5 1.882
01446500 Delaware River at Belvidere, NJ 11,745.6 1.030
01447500 Lehigh River at Stoddartsville, PA 237.5 1.633
01451000 Lehigh River at Walnutport, PA 2,302.5 0.987
01452500 Monocacy Creek at Bethlehem, PA 115.3 0.787
01463500 Delaware River at Trenton, NJ – PA side 17,560.1 1.018
01470500 Schuylkill River at Berne, PA 919.4 0.697
10 01474500 Schuylkill River at Philadelphia, PA 4,902.9 0.952
11 01467000 North Branch Rancocas Creek at Pemberton, NJ 305.6 1.031
12 01411500 Maurice River at Norma, NJ 290.1 0.975
NA NA Entire DRB (no gage) 33,288.2 1.010
No.USGS IDNameDrainage area (km2)Drainage area ratioa
01413500 East Branch Delaware River at Margaretville, NY 422.2 0.737
01415000 Tremper Kill near Andes, NY 86.0 0.601
01445500 Pequest River at Pequest, NJ 274.5 1.882
01446500 Delaware River at Belvidere, NJ 11,745.6 1.030
01447500 Lehigh River at Stoddartsville, PA 237.5 1.633
01451000 Lehigh River at Walnutport, PA 2,302.5 0.987
01452500 Monocacy Creek at Bethlehem, PA 115.3 0.787
01463500 Delaware River at Trenton, NJ – PA side 17,560.1 1.018
01470500 Schuylkill River at Berne, PA 919.4 0.697
10 01474500 Schuylkill River at Philadelphia, PA 4,902.9 0.952
11 01467000 North Branch Rancocas Creek at Pemberton, NJ 305.6 1.031
12 01411500 Maurice River at Norma, NJ 290.1 0.975
NA NA Entire DRB (no gage) 33,288.2 1.010

aRatio between the actual drainage area (shown) and the modeled drainage area (not shown).

Historic and projected future modeled land use land cover (LULC) data were obtained at 250 m resolution from the U.S. Geological Survey Land-Cover Modeling Center (Sohl et al. 2014, 2016). Data were reclassified into five LULCs: water, developed, forest, agriculture, and other. For each 12 km cell, the percentage of each of the five LULCs was calculated based on the 250 m pixels for each decade from 1950 to 2100. Interdecadal percentages for each 12 km pixel were calculated based on a linear interpolation between decadal values. Figure 1(e) shows the dominant LULC and its percentage for each cell for the year 2000. In addition to determining the amount of runoff generated by the model, percent developed is used to decrease the AWS proportionally for each grid cell.

### Model input data

The hydrologic model was developed and driven by temperature and precipitation data from the World Climate Research Programme's Coupled Model Intercomparison Project Phase 5 (CMIP5) (Taylor et al. 2012). CMIP5 data were bias-corrected and spatially downscaled to 1/8° resolution (available at: https://gdo-dcp.ucllnl.org/downscaled_cmip_projections/) (Mauer et al. 2007; Bureau of Reclamation 2013). Downscaled CMIP5 data were available on a monthly time step from January 1950 to December 2099 for 231 climate projections that were the product of 37 different climate models, 4 different RCPs (IPCC 2013), and several runs for many model/RCP combinations. See Bureau of Reclamation (2013) for a listing of all 231 projection combinations.

Also available with the downscaled climate data was a single observational dataset based on the station data of monthly temperature and precipitation from 1950 to 1999 that were gridded to the same resolution as the climate projections to allow for a comparison between the projections and observations. The observed data were used to develop the hydrologic model and the projections were used to drive the model under different conditions once the model was developed. 1950–1999 was also used as the historic baseline period to assess future change.

### Model description

The hydrologic model for the DRB is roughly based on ideas first proposed by Thornthwaite (1948). These original ideas have been expanded upon and modified in the following decades including more recently by Frei et al. (2002), Ellis et al. (2008), McCabe & Wolock (2011), Hawkins & Austin (2012), and Hawkins (2015). The model operates on a monthly time step where the processes outlined below are calculated for each grid cell in the model domain.

Monthly input precipitation totals (P) are partitioned into rainfall (Pr) and snowfall (Ps) based on values of and following comparable values in McCabe & Wolock (1999, 2011), McCabe & Markstrom (2007), and Hawkins & Austin (2012). For months when the monthly average temperature (T) is greater than Train or less than Tsnow, all precipitation is quantified as rainfall or snowfall, respectively. For months when T is between Train and Tsnow, Ps and Pr are calculated as follows:
(1)
and
(2)
Water that is generated from snowmelt (melt) and water that is stored as snow (snow) are calculated as follows:
(3)
and
(4)
where prev refers to the previous month, and smf is the snowmelt fraction and is calculated as follows:
(5)

Smf values are constrained between 0 and 1. In other studies (e.g. McCabe & Wolock 1999), smf was not allowed to exceed 0.5 to prevent an extremely large snowpack, as found in an alpine environment, from entirely melting in one warm month. For the DRB, alpine-like snowpacks do not exist and snow often entirely melts in one month. Therefore, smf was allowed to reach a value of 1.

Water available for the current month for the remainder of the hydrologic processes in the model (Pavl) is simply the sum of rainfall and snowmelt.

PE is calculated following Hamon (1961) as follows:
(6)
where d is the number of days in the month and D is the mean monthly hours of daylight in units of 12 h. If T is less than zero, PE is set to 0.

Milly & Dunne (2017) demonstrate that seven commonly used PE formulations (including the formulation in Equation (6)) that each make use of varying degrees of parameterizations in their PE estimates do a reasonable job of calculating PE when compared to the physically based (and most accurate) formulation they term ‘energy-only’. However, all seven formulations overestimate future change in PE when compared to energy-only future change estimates. Milly & Dunne (2017) suggest that this overestimate is because none of the seven methods for estimating PE account for increased plant stomatal resistance to transpiration that occurs with the increased future temperatures that climate models predict. Given that data either to calculate PE using the energy-only method or to account for the relationship between temperature and plant stomatal resistance are not available at the downscaled resolution used by the hydrologic model in this study, an empirical method to constrain the ‘runaway’, future PE calculated in Equation (6) is employed by the hydrologic model.

To constrain PE, non-downscaled near-surface air temperature and surface upward sensible (H) and latent (LE) heat flux data for the CMIP5 global climate models (GCMs) described in Milly & Dunne (2017) were obtained (https://esgf-node.llnl.gov/projects/esgf-llnl/) for the grid cell from each GCM that had the greatest intersection with the DRB. Monthly data were obtained for 1950–2099.

For each month of the data record and for each GCM, Hamon PE (PEH-GCM) was calculated as in Equation (6) and energy-only PE (PEE-GCM) was calculated following Milly & Dunne (2017) as follows:
(7)
where Rn is the net radiation at the surface, G is the heat flux into the ground, ρ is the density of water (1,000 kg/m3), and L is the latent heat of vaporization of water (2.5 × 106 J/kg). Equation (7) was converted to mm/month by multiplying by d, 86,400 s/day, and 1,000 mm/m. was evaluated by energy balance as . PEH-GCM and PEE-GCM values were averaged by month, RCP, and for different time slices over all the GCMs.
To constrain the hydrology model PE, PE for each month and for each hydrology model grid cell was compared to the historical average of PEH-GCM for the appropriate month and RCP. If the difference between the hydrology model PE and the historical average of PEH-GCM was greater than that between PEE-GCM for the appropriate month, year, and the RCP and the historical average of PEE-GCM, then PE was constrained and recalculated as follows:
(8)
where and are the 1950–1999 average PE and values, respectively, and is the 10-year running mean of monthly PE. This recalculation of PE constrains ‘runaway’, temperature-based PE based on the results in Milly & Dunne (2017), while at the same time preserving the time-series variability of PE calculated in Equation (6).
Subsurface moisture storage (SM) is calculated as follows:
(9)
where . If SM exceeds the AWS, surplus water (surp) is generated as the difference between SM and AWS. If Pavl ≥ PE, then actual evapotranspiration (AE) equals PE. Otherwise, AE equals Pavl minus the change in SM from the previous month. Runoff (RO) is calculated as follows:
(10)
where R is the runoff ratio that determines how much surplus runs off the current month and how much is held over until the next month. R is a weighted average value based on the percentage LULCs for each grid cell and the R-values for those LULCs show in Table 2. Note that because R is a function of LULC which changes with time, R changes slightly with time. Grid cell runoff is also aggregated for the subwatersheds in Figure 1(d).
Table 2

Runoff coefficients (R) for different LULCs used in the hydrology model

LULCR
Water 0.9
Developed 0.8
Forest 0.5
Agriculture 0.5
Other 0.7
LULCR
Water 0.9
Developed 0.8
Forest 0.5
Agriculture 0.5
Other 0.7

Values are comparable to those reported previously (Mather 1978; McCabe & Wolock 2011; Hawkins & Austin 2012).

### Model evaluation

Monthly and water year annual runoff from the hydrologic grid cells were aggregated to the subwatersheds shown in Figure 1(d) and Table 1. Modeled aggregated subwatershed runoff was multiplied by the ratio between the subwatershed drainage area and the area comprised by the model grid cells that represented the subwatershed. Due to the coarse subwatershed boundaries associated with the 12 km grid cells, modeled and actual watershed areas often differed, especially for the smaller subwatersheds. Modeled runoff values were compared with observed runoff obtained from the USGS National Water Information System: Web Interface (https://waterdata.usgs.gov/nwis). Standard evaluation statistics were calculated that included: Pearson correlation (r), root-mean-squared error (RMSE), mean absolute error (MAE), and Nash–Sutcliffe (NS) efficiency coefficient (Nash & Sutcliffe 1970).

## RESULTS

### Model evaluation

Figure 2 shows water year annual and monthly model-predicted and -observed stream flow for the two largest rivers in the DRB, the Delaware and Schuylkill Rivers, as well as predicted runoff for the entire DRB. Similar graphs for the other gages in Table 1 are not shown. Table 3 shows model evaluation statistics for all subwatersheds. Graphically, water year predictions match observations nearly perfectly. Quality is further demonstrated by correlation coefficients that range from 0.884 to 0.981 with an average value of 0.944 and NS efficiency coefficients that range from 0.416 to 0.962 with an average value of 0.721. An NS value of 1 corresponds to a perfect match between observations and predictions, while a value of 0 indicates that the predictions were as accurate as using the mean of the observations. A negative value indicates that the mean of the observations was a better predictor. The model performed better on the larger watersheds. MAE and RMSE were normalized by dividing by the mean flow to allow comparisons between watersheds where the flow varies significantly. Again, based on the error statistics, the larger watersheds had lower errors. Modeled water year average runoff for the entire DRB ranged from 236 to 833 m3/s and averaged 560 m3/s between 1950 and 1999.

Table 3

Model evaluation statistics for comparing predicted and observed runoff for the 12 subwatersheds on the DRB on a monthly and water year timescale.

Monthly
Water year
Gage no.rMAE/meanRMSE/meanNSrMAE/meanRMSE/meanNS
0.807 0.40 0.54 0.637 0.950 0.09 0.12 0.796
0.805 0.36 0.52 0.621 0.894 0.11 0.13 0.658
0.870 0.33 0.44 0.628 0.967 0.19 0.22 0.559
0.882 0.25 0.35 0.777 0.971 0.06 0.07 0.934
0.856 0.27 0.41 0.707 0.945 0.12 0.13 0.665
0.881 0.25 0.36 0.742 0.956 0.13 0.14 0.711
0.770 0.56 0.68 −0.365 0.891 0.22 0.28 0.416
0.896 0.22 0.31 0.803 0.981 0.04 0.05 0.962
0.895 0.25 0.37 0.759 0.977 0.12 0.14 0.768
10 0.912 0.24 0.32 0.824 0.972 0.08 0.09 0.918
11 0.855 0.39 0.47 0.204 0.939 0.10 0.12 0.816
12 0.831 0.42 0.50 −0.056 0.884 0.17 0.21 0.447
Average 0.855 0.33 0.44 0.523 0.944 0.12 0.14 0.721
Monthly
Water year
Gage no.rMAE/meanRMSE/meanNSrMAE/meanRMSE/meanNS
0.807 0.40 0.54 0.637 0.950 0.09 0.12 0.796
0.805 0.36 0.52 0.621 0.894 0.11 0.13 0.658
0.870 0.33 0.44 0.628 0.967 0.19 0.22 0.559
0.882 0.25 0.35 0.777 0.971 0.06 0.07 0.934
0.856 0.27 0.41 0.707 0.945 0.12 0.13 0.665
0.881 0.25 0.36 0.742 0.956 0.13 0.14 0.711
0.770 0.56 0.68 −0.365 0.891 0.22 0.28 0.416
0.896 0.22 0.31 0.803 0.981 0.04 0.05 0.962
0.895 0.25 0.37 0.759 0.977 0.12 0.14 0.768
10 0.912 0.24 0.32 0.824 0.972 0.08 0.09 0.918
11 0.855 0.39 0.47 0.204 0.939 0.10 0.12 0.816
12 0.831 0.42 0.50 −0.056 0.884 0.17 0.21 0.447
Average 0.855 0.33 0.44 0.523 0.944 0.12 0.14 0.721

MAE and RMSE are normalized dimensionless values that are calculated by dividing the MAE and RMSE by the mean observed flow.

Figure 2

Water year annual (a and b) and monthly (d and e) observed and predicted stream flow for gage 8, Delaware River at Trenton, NJ – Pennsylvania side, and gage 10, Schuylkill River at Philadelphia, PA (Figure 1 and Table 1). Also, predicted streamflow for the entire DRB (c and f). All data are 1950–1999.

Figure 2

Water year annual (a and b) and monthly (d and e) observed and predicted stream flow for gage 8, Delaware River at Trenton, NJ – Pennsylvania side, and gage 10, Schuylkill River at Philadelphia, PA (Figure 1 and Table 1). Also, predicted streamflow for the entire DRB (c and f). All data are 1950–1999.

Monthly average predictions do not graphically align with observations as well as water year values. While monthly statistics are not as strong as annual statistics, the model still performs well on a monthly basis. The model captures the timing of the annual peak and minimum discharge. The model matches observed values better during the warm season from May to October compared with the cold season from November to April. Correlations range from 0.770 to 0.912 with an average value of 0.855. NS values ranged from −0.365 to 0.824 with an average value of 0.523. MAE and RMSE values for every watershed are larger for the monthly average assessment compared to the water year assessment. Again, the model performed better on larger watersheds. For the entire DRB, streamflow on average from 1950 to 1999 reached a maximum in March at 967 m3/s and a minimum of 196 m3/s in August.

### Basinwide model results

Figure 3 shows water year and monthly average hydrologic output averaged by the RCP and time slices out to 2099 for the entire DRB. Figure 4 shows the percent change between hydrologic outputs (excluding the temperature) between future time slices and the baseline 1950–1999 period, and Table 4 shows the trends per century for these same data. The annual average temperature is projected to increase 1.9–4.5 °C/century, depending on the RCP (Figure 3(a) and Table 4). Increases are relatively constant for each month with slightly greater warming occurring during winter and summer (Figure 3(b)).

Table 4

Trends per century for hydrologic variables for the entire DRB averaged by the RCP

RCP 2.6RCP 4.5RCP 6.0RCP 8.5
T (°C) 1.9 2.8 3.1 4.5
PE (mm) 58.4 78.4 107.1 99.5
AE (mm) 49.1 64.5 84.0 76.7
SM (%) −2.1 −3.3 −5.4 −5.0
RO (mm) 48.9 51.5 41.0 72.7
P (mm) 98.1 116.3 125.6 150.0
Pr (mm) 132.4 162.1 172.1 212.6
Ps (mm) −34.3 −45.7 −46.5 −62.5
Sn (mm) −3.9 −4.8 −4.8 −5.9
M (mm) −34.3 −45.7 −46.5 −62.5
RCP 2.6RCP 4.5RCP 6.0RCP 8.5
T (°C) 1.9 2.8 3.1 4.5
PE (mm) 58.4 78.4 107.1 99.5
AE (mm) 49.1 64.5 84.0 76.7
SM (%) −2.1 −3.3 −5.4 −5.0
RO (mm) 48.9 51.5 41.0 72.7
P (mm) 98.1 116.3 125.6 150.0
Pr (mm) 132.4 162.1 172.1 212.6
Ps (mm) −34.3 −45.7 −46.5 −62.5
Sn (mm) −3.9 −4.8 −4.8 −5.9
M (mm) −34.3 −45.7 −46.5 −62.5

T, temperature; PE, potential evapotranspiration; AE, actual evapotranspiration; SM%, subsurface moisture percentage; RO, runoff; P, precipitation; Pr, rainfall; Ps, snowfall; Sn, snow water storage; M, snowmelt.

Figure 3

Entire DRB water year annual average or total (columns 1 and 3) and monthly average (columns 2 and 4) data from the hydrologic model when driven with future climate projections. Water year data are averaged by the RCP. Monthly data are averaged by the RCP and time slices.

Figure 3

Entire DRB water year annual average or total (columns 1 and 3) and monthly average (columns 2 and 4) data from the hydrologic model when driven with future climate projections. Water year data are averaged by the RCP. Monthly data are averaged by the RCP and time slices.

Figure 4

Percent change in hydrologic variables for the entire DRB, averaged by RCP, between 2020–2039, 2050–2069, and 2080–2099 compared to the 1950–1999 baseline. PE, potential evapotranspiration; AE, actual evapotranspiration; SM, subsurface moisture; RO, runoff; P, precipitation; Pr, rainfall; Ps, snowfall; Sn, snow water storage; M, snowmelt.

Figure 4

Percent change in hydrologic variables for the entire DRB, averaged by RCP, between 2020–2039, 2050–2069, and 2080–2099 compared to the 1950–1999 baseline. PE, potential evapotranspiration; AE, actual evapotranspiration; SM, subsurface moisture; RO, runoff; P, precipitation; Pr, rainfall; Ps, snowfall; Sn, snow water storage; M, snowmelt.

Annual average precipitation is projected to increase 98.1–150.0 mm/century, depending on the RCP (Figure 3(c) and Table 4). This equates to an 8–16% increase by 2080–2099 (Figure 4). Increases are relatively constant across each month with slightly greater increases in precipitation occurring during winter and spring (Figure 3(d)). Due to warming, the fraction of precipitation that falls as rainfall is projected to increase 132.4–212.6 mm/century (Table 4) or 12–24% by 2080–2099 (Figures 3(g) and 4), while snowfall is projected to decrease 34.3–62.5 mm/century (Table 4) or 42–84% by 2080–2099 (Figures 3(g) and 4). Unsurprisingly, the biggest changes in rainfall and snowfall are projected to occur in winter, with a maximum change in January, historically the snowiest month, as warmer temperatures convert snowfall to rainfall (Figure 3(l)). Under RCP 8.5 by 2080–2099, basinwide DJF snowfall is projected to decrease 13.1 mm, an 84% decrease from the 80.0 mm 1950–1999 baseline.

Correspondingly, snow water storage is projected to decrease 3.9–5.9 mm/century (Table 4), a 58–93% decrease by 2080–2099 (Figures 3(o) and 4), and snowmelt is projected to decrease 34.3–62.5 mm/century (Table 4), a 42–84% decrease by 2080–2099 (Figures 3(s) and 4). Again, changes in snow water storage and snowmelt are projected to be greatest in winter (Figures 3(p) and 3(t)), with the largest change in snow water storage occurring in January and February, historically the months of deepest snowpack. Interestingly, the snowmelt peak is projected not only to decrease but also to shift earlier in the season. Historically, snowmelt peaked in March. This remains true for RCP 2.6 in 2020–2039. However, for RCP 8.5 in 2020–2039, both RCPs in 2050–2069 and RCP 2.6 in 2080–2099, the snowmelt peak shifts to February. For RCP 8.5 in 2080–2099, the snowmelt peak shifts to January.

Due to the projected warming, PE is projected to increase 58.4–107.1 mm/century (Table 4) or 9–19% by 2080–2099 (Figures 3(e) and 4), while AE is projected to increase 49.1–84.0 mm/century (Table 4) or 8–16% by 2080–2099 (Figures 3(i) and 4). Projected increases in PE and AE on a monthly basis are greatest in the warmer summer months and minimal during the colder winter months (Figures 3(f) and 3(j)). During the 1950–1999 baseline period, PE demand was not met by AE from May through October. This resulted in utilization and conversion to AE of subsurface moisture to attempt to meet the demand. This timeframe has mostly been projected to remain the same with a few RCP/time slice combinations shifting the deficit period to April through October. Historically, the annual average PE–AE deficit was 20.8 mm. This deficit increases to 29.6–47.2 mm (42–128% increase) for different RCP/time slice combinations.

As a result of this increasing PE–AE deficit, SM is projected to decrease 2.1–5.0%/century (Table 4) or 3–7% by 2080–2099 (Figures 3(m) and 4). Note that SM is presented as % capacity. Therefore, percentage values are interpreted differently than percent change for data measured in units of degrees or mm. During the 1950–1999 baseline period, SM was below 99% capacity from May through November (Figure 3(n)). Subsurface moisture was recharged above 99% by December. Note that values are basinwide averages and thus certain cells may not fully recharge. Therefore, 99% was chosen to represent full recharge when averaged over the entire DRB. This timeframe has mostly been projected to remain the same with a few RCP/time slice combinations shifting the beginning of the SM utilization period to April. Historically, percent SM capacity is the lowest in July at 63.6%. For RCP 8.5 in 2080–2099, percent capacity drops to 52%, an 18% decrease.

Finally, runoff, as an integrator of the previously discussed variables, is projected to increase 41.0–72.7 mm/century (Table 4) or 8–17% by 2080–2099 (Figures 3(q) and 4). Monthly runoff is essentially unchanged from March through November (Figure 3(r)). During the winter months, runoff is projected to increase as more precipitation falls as rain instead of snow, and snowmelt increases and occurs earlier than during the baseline period. Winter (DJF) total increases in runoff depth range from 25.6 to 74.0 mm, an increase of winter runoff of 15–43%.

### Gridded model results

Figures 5 and 6 show the geographic pattern of changes for the same variables that are presented in Figure 3 as watershed average values. For temperature (Figure 5(a)), all grid cells for all cases show warming. The pattern of warming is consistent across all maps with greater warming at more inland locations that are farther removed from the moderating influence of the Atlantic Ocean. These inland locations are at higher elevation and generally in the north. The range in differences between the most and least warming for each map is from 0.1 to 1.4 °C. The degree of warming is fairly consistent across seasons. However, spring showed slightly less warming in all cases. The greatest warming occurs in the mid- and late century under RCP 8.5. Basinwide average temperature increases for RCP 8.5 by 2080–2099 are greatest in summer at 5.8 °C and least in spring at 4.9 °C.

Figure 5

Gridded hydrologic model data maps. Rows separated by heavy dark lines are (a) temperature, (b) potential evapotranspiration, (c) actual evapotranspiration, (d) percent capacity of subsurface moisture, and (e) runoff. For each variable, the top row represents RCP 2.6 and the bottom row represents RCP 8.5. For each variable, the columns represent four seasonal maps of changes between 2020–2039, 2050–2069, and 2080–2099 when compared with the baseline 1950–1999. Because the range of values differed greatly between RCPs and time slices, the color legend is a relative scale and colored values within the watershed boundaries are different for each map. This allows for the geographic variability to be seen but does not allow for comparison between maps. To compare maps, the black and white legend represents an absolute scale of the basinwide average of the gridded data. This legend applies to the background shading outside the watershed boundary for each map. The average value of the gridded data on which the black and white shading is based is presented in the top left corner of each map. For each set of four seasonal maps, the largest average value is bold and the smallest value is italicized. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wcc.2019.140.

Figure 5

Gridded hydrologic model data maps. Rows separated by heavy dark lines are (a) temperature, (b) potential evapotranspiration, (c) actual evapotranspiration, (d) percent capacity of subsurface moisture, and (e) runoff. For each variable, the top row represents RCP 2.6 and the bottom row represents RCP 8.5. For each variable, the columns represent four seasonal maps of changes between 2020–2039, 2050–2069, and 2080–2099 when compared with the baseline 1950–1999. Because the range of values differed greatly between RCPs and time slices, the color legend is a relative scale and colored values within the watershed boundaries are different for each map. This allows for the geographic variability to be seen but does not allow for comparison between maps. To compare maps, the black and white legend represents an absolute scale of the basinwide average of the gridded data. This legend applies to the background shading outside the watershed boundary for each map. The average value of the gridded data on which the black and white shading is based is presented in the top left corner of each map. For each set of four seasonal maps, the largest average value is bold and the smallest value is italicized. Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/wcc.2019.140.

Figure 6

Same as Figure 5 except rows separated by heavy dark lines are (a) precipitation, (b) rainfall, (c) snowfall, (d) snow water storage, and (e) snowmelt.

Figure 6

Same as Figure 5 except rows separated by heavy dark lines are (a) precipitation, (b) rainfall, (c) snowfall, (d) snow water storage, and (e) snowmelt.

The pattern of the precipitation change is not as consistent as temperature change (Figure 6(a)). However, seasonal patterns are fairly consistent across RCPs and time slices. All cells for all cases show increases in precipitation. The range in differences between the smallest and largest increases in precipitation for each map is from 1.9 to 11.8 mm. Winter shows the greatest increases in the upper three quarters of the DRB with the exception of the northernmost tip of the watershed where increases are slightly less. Spring shows a similar pattern but the northernmost tip does not register smaller increases in precipitation in the manner that winter did. With one exception, summer precipitation increases are greatest in the southern part of the watershed and begin a gradient toward smaller increases in the north starting near the midpoint of the north–south transect. RCP 2.6 for 2020–2039 in summer does not have this north–south gradient but rather has smaller increases in the westernmost part of the watershed. Fall shows fairly uniform precipitation increases across the southern two-thirds of the watershed. Early in the century, fall shows the greatest increases in the northernmost part of the watershed and smallest increases immediately south of the largest increase area. By the end of the century, the greatest fall increases are only seen in the most northeastern part of the watershed as the rest of the upper third is showing smaller increases. Average watershed increases are generally greater for RCP 8.5 and later in the century. For all RCP 8.5 time slices and for RCP 2.6 in 2020–2039, the greatest average increases in precipitation occur in winter. For all RCPs and time slices, the smallest average increases occur in fall.

Figures 6(b) and 6(c) elaborate on the precipitation changes by showing patterns of changes in rainfall and snowfall. Summer rainfall and precipitation maps are identical and summer snowfall maps have no data as all historical summer precipitation is rainfall. All grid cells for all cases show increases in rainfall and decreases in snowfall. The range in differences between the smallest and largest increases in rainfall for each map is from 2.2 to 37.3 mm, while single map ranges in snowfall decreases are from 2.2 to 47.0 mm. Ranges are the greatest in winter and the smallest in fall. In all cases, rainfall increases and snowfall decreases are greatest in winter as the warmer temperatures convert precipitation from snowfall to rainfall. These changes are most pronounced in the upper two-thirds of the watershed where the greatest amount of snow has fallen historically. At the most northern and high elevation area of the DRB in winter, increases in rainfall and decreases in snowfall are slightly less than the rest of the upper basin as these areas also had slightly smaller increases in overall precipitation (Figure 6(a)) in winter. Patterns of spring and fall rainfall increases are similar to precipitation increases, but there is a greater increase in rainfall compared to precipitation for the upper quarter of the DRB.

Correspondingly, snowfall decreases for spring and fall are greatest in the upper quarter of the DRB as these are the areas that have traditionally had snowfall that has been converted to rainfall during the winter shoulder seasons (mostly March and to a lesser extent April and November). Fall basinwide average rainfall increases and snowfall decreases are smallest (excluding summer). Of the three snowfall seasons, fall has the lowest historic snowfall. Overall, changes in rainfall and snowfall are greatest for RCP 8.5 and later in the century.

Patterns of changes in snow water storage and snowmelt (Figures 6(d) and 6(e)) closely mirror patterns of changes in snowfall. Again, there are no summer data. Snow water storage for all cells for all cases shows decreases. The range of values in the fall maps is close to zero as all changes are less than 1 mm as snow typically melts in the fall and is not stored. Winter and spring ranges in snow water storage decreases are from 33.1 to 93.0 and 14.8 to 25.9 mm, respectively, and all ranges increase moving toward the end of the century. Decreases in snow water storage are greatest in the northern and high elevation areas of the watershed where snowpack has traditionally existed and where warming is greatest. Basinwide decreases in snow water storage are greatest in winter, much less in spring, and nearly zero in fall. Figure 7 shows the projected snow water storage that results from the changes in snowfall, snow water storage, and snowmelt. A dramatic decrease can be seen for both RCPs with snow water storage nearly disappearing under RCP 8.5 by 2080–2099.

Figure 7

DJF average snow water storage for one historic and three future periods and for RCPs 2.6 and 8.5. Cells with values less than 10 mm are not shown. For each map, shown in the upper left corner are the number of cells with snow water storage greater than 10 mm and the average snow water storage of those cells.

Figure 7

DJF average snow water storage for one historic and three future periods and for RCPs 2.6 and 8.5. Cells with values less than 10 mm are not shown. For each map, shown in the upper left corner are the number of cells with snow water storage greater than 10 mm and the average snow water storage of those cells.

Snowmelt for all cells in spring and fall for all cases shows decreases with the greatest decreases in the north and high elevations. For winter, snowmelt increases in the northernmost part of the DRB and decreases for the rest of the DRB. The area of increased snowmelt decreases as the century progresses and disappears under RCP 8.5 by the end of the century. Increased winter snowmelt corresponds to the combination of warming temperatures but still enough snowpack available to melt. Spring, winter, and fall ranges in snowmelt changes are from 21.2 to 57.2, 13.9 to 19.8, and 1.6 to 2.4 mm, respectively, and all ranges increase moving toward the end of the century.

For PE (Figure 5(b)), all grid cells for all cases show increases. The range of values for every map is less than 1 mm, indicating that while these changes can be mapped, there are no meaningful differences in PE changes across the basin for each map. While there is a clear pattern to temperature changes across the basin (Figure 5(a)), the influence of greater warming in the north is likely offset by the increasing latitude that factors into the PE calculation as well as the impact of constraining the PE based on GCM grid cells with differing geographies. For both RCPs presented, increases in PE are greatest in summer and smallest in winter with a summer 2080–2099 RCP 8.5 basinwide average PE increase of 15.1 mm.

Patterns of changes in AE (Figure 5(c)) closely align with PE patterns. In fact, for winter for all cases, the PE and AE maps are identical (presented in Figure 5 with different colors). With the exception of summer, AE increases for all cells for all cases. For summer, nearly all of the cells show increased AE. Summer ranges in AE changes are from 16.2 to 23.1 mm. Fall and spring ranges are from 0.9 to 4.0 mm, and winter ranges are near zero. All ranges increase moving toward the end of the century. Summer AE decreases and spring and fall smaller increases occur in the same geographic area that corresponds to smaller natural subsurface moisture storage capacity (Figure 1(c)). This smaller capacity is further decreased due to the urbanization associated with Philadelphia. The demand for water due to temperature-driven increases in PE cannot be met with precipitation and subsurface moisture in areas where the subsurface storage is smaller, thus resulting in decreased AE in summer and smaller increases in AE during spring and fall. For all RCP 2.6 time slices and for RCP 8.5 in 2020–2039, basinwide average AE shows the greatest increase in summer due to the increased PE and the ability of the subsurface to still supply water to AE. For mid- and late-century RCP 8.5, basinwide AE increases are greatest in spring due to the availability of subsurface moisture which has been diminished by summer.

Figure 5(d) shows changes in subsurface moisture storage. The majority of grid cells show decreases in SM. The notable exception is winter where the northern part of the watershed shows increases in SM. In winter, with low PE and ample precipitation, nearly every cell for every year is at capacity, so changes in SM are not particularly meaningful. The range in winter SM changes is about 1 mm. Small winter increases are the result of any unfilled capacity being filled by increased rainfall and snowmelt. The southern part of the watershed is characterized by decreases in SM that are characteristic of warming temperatures and increased PE. Summer ranges in SM changes are from 25.8 to 33.4 mm. Fall and spring ranges are from 7.1 to 13.0 mm. Basinwide average decreases are greatest in summer when PE demand is the highest and SM is utilized for AE. Decreases are the second largest in fall. Changes are most pronounced in the urban Philadelphia area where storage capacity is also naturally the smallest. Otherwise, SM decreases across the DRB for spring, summer, and fall are fairly uniform for each individual map. The very few spring, summer, and fall grid cells that show slight increases in SM occur in urban areas with low natural capacity. Therefore, a few wet years are able to quickly recharge these cells and pull the average toward increases.

Figure 8 shows the JJA SM values that result from the changes. Summer is the minimum for SM. Summer SM decreases with time and is least under RCP 8.5. Lowest SM values are found near urban areas with small natural storage capacity. Utilization of SM due to AE demand during the warm summer months is most impactful on the small capacity cells.

Figure 8

JJA average subsurface moisture percentage for one historic and three future periods and for RCPs 2.6 and 8.5. For each map, shown in the upper left corner are the number of cells with subsurface moisture percentage less than 50% and the average subsurface moisture percentage of all the cells on the DRB.

Figure 8

JJA average subsurface moisture percentage for one historic and three future periods and for RCPs 2.6 and 8.5. For each map, shown in the upper left corner are the number of cells with subsurface moisture percentage less than 50% and the average subsurface moisture percentage of all the cells on the DRB.

Finally, Figure 5(e) shows patterns in runoff change. Basinwide average runoff shows increases for all cases. Winter is the only season when runoff increases for all cells in all cases. Winter is also the season with the greatest increases, with ranges in runoff changes from 17.7 to 43.1 mm. The greatest increases occur in the northern high elevation areas where snowfall is converting to rainfall, snow water storage is decreasing, and snowmelt is increasing. The other three seasons show a general pattern of decreased runoff in the north and increased runoff in the south. This pattern is most pronounced in spring, and to a lesser extent fall, due to decreased snow water storage and snowmelt as these phenomena increasingly occur earlier in the year. Summer decreases in runoff are due to increased AE as well as smaller increases in rainfall in the north compared to the rest of the DRB. Ranges for spring and summer runoff differences are from 5.2 to 19.0 mm. Ranges for fall differences are from 15.5 to 22.7 mm. All ranges are larger with time and for RCP 8.5.

## DISCUSSION

### Practical implications

As a drinking, agricultural, and industrial water supply to over 15 million people, changes in the timing and amount of runoff for the DRB have significant implications. Currently, DRB reservoir operations are governed by the Flexible Flow Management Program (USGS 2017), the Delaware River Basin Water Code (DRBC 2013), and the DRBC Drought Operating Plans (DRBC 1983, 1988). Changes in reservoir operations and water allocations may need to be considered as winter runoff increases due to increased precipitation, increased winter rainfall:snowfall ratio, increased snowmelt, and low winter evaporation. Historically, DRB runoff has peaked in March but in the future could peak in January or February depending on the geographic location, RCP, and future time slice. Consideration of reservoir water supply storage in the fall will become of greater concern if peak runoff occurs months earlier than it has historically. Also, as future flooding is considered (Woltemade 2016; Woltemade et al. in review), reservoir management as related to flood control should be assessed.

While the above changes in the timing of water supply are critical to plan for, it is equally important to consider changes in water demand as the DRB warms. Traditionally, evaporation from water supply reservoirs has been a larger concern for reservoirs in the western USA (e.g. Friedrichi et al. 2018). However, as eastern US watersheds, including the DRB, warm it may be worth considering the impact that increased temperatures and therefore increased PE have on water supplies. This is especially true as traditionally cooler months begin to see increased evaporation and inflow into reservoirs occurs earlier.

Additionally, on the demand side of the DRB water resource equation, decreasing subsurface moisture during the summer under a warming climate may become a larger concern, especially as related to agriculture in the lower basin and forest health in the upper basin. This concern is especially problematic during drought conditions when subsurface moisture deficits could be significantly enhanced. The DRBC has declared six drought emergencies since the 1960s (DRBC 2019b). The impact of a warmer climate on future droughts deserves serious consideration so as to be adequately prepared.

### Hydrologic model uncertainty

#### Input temperature and precipitation data

As noted by many previous studies of DRB hydrology, hydrologic models of future conditions are dependent on the quality of the input climate data used to drive the model. As with all of the hydrologic studies referenced in this report, the variability associated with projections of precipitation for this study is much larger than temperature as the physics associated with modeling precipitation often occurs at a resolution smaller than that of most climate model grid cells. As precipitation is the main driver runoff in the DRB, uncertainty in precipitation input data affects runoff projections. This uncertainty is further compounded by the fact that for the DRB, future runoff changes are also highly dependent on changes in snowfall, snowmelt, and snow water storage, all of which are functions of input precipitation data. While uncertainty in precipitation data is greater than temperature, results from this study are consistent with what has been reported (e.g. Najjar et al. 2009) for years for the mid-Atlantic region: annual precipitation is projected to increase with the greatest increases occurring in winter. Thus, there is confidence that the ensemble precipitation projection is a meaningful hydrologic variable despite the larger range in the individual precipitation projections compared to temperature.

#### Model output distribution

Figure 9 shows the distribution of percent change in runoff for all 231 individual climate projections. As previously discussed in Figure 3, basinwide runoff increases with larger forcing RCPs and at times further in the future. The standard deviation associated with the distributions also generally increases with larger forcing RCPs and absolutely increases at times further in the future, suggesting that the larger changes in climate associated with larger forcing RCPs and times further in the future have a compounding effect on the output of the hydrologic model.

Figure 9

Frequency distribution of % difference in the DRB aggregate runoff by the RCP for (a) 2020–2039, (b) 2050–2069, and (c) 2080–2099 compared to the 1950–1999 baseline for 231 climate projections used to drive the hydrologic model. Vertical dotted lines indicate the RCP average value. Also, (d) shows the standard deviation of the distributions by the RCP and time slice.

Figure 9

Frequency distribution of % difference in the DRB aggregate runoff by the RCP for (a) 2020–2039, (b) 2050–2069, and (c) 2080–2099 compared to the 1950–1999 baseline for 231 climate projections used to drive the hydrologic model. Vertical dotted lines indicate the RCP average value. Also, (d) shows the standard deviation of the distributions by the RCP and time slice.

#### Constrained versus unconstrained PE

Figure 10 shows some of the results of running the hydrologic model using constrained versus unconstrained PE. By 2100, constrained PE is 21–180 mm/year less than unconstrained PE depending on the RCP (Figures 10(a) and 10(b)). These differences are greatest in summer (not shown). As discussed in the Methods section, constrained PE was chosen for the hydrologic model. The resulting AE by 2100 was 13–105 mm less than AE based on unconstrained PE. Subsurface moisture percent capacity was 2–11% greater and runoff was 14–104 mm greater for constrained PE compared to unconstrained PE. The impact of constrained versus unconstrained PE on runoff is relatively minimal for RCPs 2.6, 4.5, and 6.0 as runoff increases slightly using constrained PE (Figure 10(d)), but the nature of the runoff over time (increasing, decreasing, or flat) is not impacted (not shown). For RCP 8.5, not only is runoff greater using constrained PE (Figure 10(d)) but the nature of runoff with time changes (not shown). Constrained PE produces increased runoff by 2100 compared to historical baselines, while unconstrained PE produces slightly decreasing runoff. As discussed in Milly & Dunne (2017), PE estimates, such as the unconstrained values in this study, overestimate future PE and therefore can have profound implications on other components of the hydrologic cycle as illustrated here. It is thus critical to evaluate the impacts of PE on hydrologic models and employ a methodology to limit ‘runaway’, unconstrained PE as was done in this study. Williamson et al. (2016) found similar results for the DRB in regards to PE formulation and runoff.

Figure 10

Constrained and unconstrained water year annual PE totals (a). Also, differences in (b) PE and AE, (c) subsurface moisture, and (d) runoff that resulted from the hydrologic model being run with constrained or unconstrained PE.

Figure 10

Constrained and unconstrained water year annual PE totals (a). Also, differences in (b) PE and AE, (c) subsurface moisture, and (d) runoff that resulted from the hydrologic model being run with constrained or unconstrained PE.

### Model performance

#### Subwatershed size

When comparing observed and predicted streamflow, the hydrologic model performed better on the larger subwatersheds that were evaluated. While a ratio between the actual and modeled watershed area was used to correct the modeled, aggregated streamflow values, ratios strayed farther from 1.0 for the smaller watersheds. Additionally, for the smaller watersheds, in some cases, less than four grid cells were aggregated to produce a watershed streamflow value. This does not imply that calculated values from any given grid cell are more accurate than another grid cell. Rather, the scale and geography of the smaller watersheds do not align as well with the observed data and therefore make comparison less meaningful. The model performed very well on the larger watersheds which are the focus of this study.

#### Timescale

The model also performed better when examining water year values compared to monthly average values. Water year stream flow is highly dependent on the input of precipitation and slightly less dependent on the outflow of water due to evapotranspiration. Other processes within the model serve to change the timing of the monthly hydrologic processes, but, at the water year resolution, these processes do not matter when compared to precipitation and evapotranspiration. Given that precipitation changes on a daily or even sub-daily timescale, it could be expected that the monthly resolution predictions would perform as well as the annual predictions. This would likely be true in a warmer climate or a snow-dominated winter climate. For the DRB however, the sub-monthly resolution of actual snowfall, snow accumulation, and snowmelt is not captured well by the monthly time step of the model. In the DRB, it happens regularly that snow falls, accumulates, and melts all within the same month or perhaps multiple times. More problematic from the perspective of a model with a monthly time step is snow that falls, accumulates, or quickly melts, but these processes occur in different months. Similarly, problematic is the fact that the model uses a monthly average temperature to determine precipitation type, snow accumulation, and snowmelt. It is likely that there are months when the average monthly temperature dictates a certain type of precipitation and the amount of snowmelt, but actual individual precipitation or melt events have associated temperatures that deviate significantly from the monthly average temperature. All of these sub-monthly processes serve to skew the runoff generated by a model with a monthly time step. It is likely that monthly runoff values would match better with observations in a warmer watershed where snowfall does not occur or a higher elevation, mountainous watershed where snow accumulates throughout fall and winter and has a steady, consistent melt in the spring.

Given these monthly limitations, the model still does a satisfactory job of capturing the timing of the peak and low flow of the annual runoff cycle as well as the overall shape of the annual hydrograph. Similar to the water year assessment and for the same reasons, performance is better for the larger watersheds. For the entire DRB, and for most of the subwatersheds examined, runoff increases from September through December. This is primarily due to a decrease in evapotranspiration as the weather cools and a nearly constant precipitation amount that falls as rain. From December to February, runoff stays nearly constant as much of the precipitation now falls as snow and the water is stored on the surface as snowpack rather than generating runoff. From February to April, runoff increases again due to a combination of increased snowmelt and precipitation that is falling more as rain. Finally, for the remainder of the year, runoff decreases as evapotranspiration increases with the warmer months, and more subsurface storage is available to absorb potential runoff.

## CONCLUSION

A gridded hydrological model was developed for the DRB that was driven by future climate projections to assess the impact on DRB runoff and other hydrologic conditions. A novel method for calculating PE is presented. Major findings of this study are summarized as follows:

• Basinwide, annual temperature, PE, AE, precipitation, rainfall, and runoff are all projected to increase while snowfall, snow water storage, snowmelt, and subsurface moisture are all projected to decrease. Temperature increases driving these changes range from 2.0 to 5.5 °C by 2080–2099.

• As integrators of all the hydrologic processes, basinwide summer subsurface moisture is projected to decrease 7–18% due to increased PE, while winter runoff is projected to increase 15–43% due to increased precipitation/rainfall and increased snowmelt.

• Significant spatial variability in future changes to all hydrologic conditions exists across the DRB. Of particular note, winter changes in rainfall, snowmelt, and runoff are most pronounced in the traditionally snow covered regions of the high elevation, northern DRB.

While the results of this study are specific to the DRB, several key points have broader general applicability: (1) substantial spatial variability exists in the hydrologic response to climate change across regional river basins of the scale investigated here, (2) resource managers should expect to plan for future water resource conditions that will vary from the past due to hydroclimatologic changes, (3) a few key hydrologic variables are likely to play dominant roles in the impact of future climate conditions on water resources, including soil moisture, PE, and snowmelt, and (4) uncertainty in climate projections warrants added flexibility and adaptive management approaches to future water resource management.

Ongoing and future work will address the implications of future climate and hydrologic changes on flooding, drought, and forest ecosystem health for the DRB. Additionally, the hydrologic model will be run using higher resolution climate input data to assess model quality as a function of resolution. Collectively, these projects will aid water resource managers as they plan for future infrastructure and operations procedures to manage the watershed.

## ACKNOWLEDGEMENTS

This work was funded by grant #DWRF-16-106 from the Delaware Watershed Research Initiative. We thank to Alfonso Yáñez Morillo, Jonathon Chester, and Troy Alleman for help processing GIS data and also thank to Claire Jantz, Scott Drzyzga, Patrick Jantz, and Antonia Price for insightful comments.

## REFERENCES

REFERENCES
Andrus
A.
Roberts
E.
Dalke
S. P.
(eds)
2017
.
Pinchot Institute for Conservation
,
Washington, DC
.
Bastola
S.
2013
Hydrologic impacts of future climate change on the southeast US watersheds
.
Regional Environmental Change
13
(
Suppl. 1
),
131
139
.
https://doi.org/10.1007/s10113-013-0454-2
.
Beecher
S.
Thaler
T.
Griffith
G.
Perry
A.
Crossett
T.
R.
2013
Adapting to A Changing Climate: Risks and Opportunities for the Upper Delaware River Region
.
Model Forest Policy Program in Association with Common Waters Partnership, Pinchot Institute for Conservation, The Cumberland River Compact, and Headwaters Economics
,
Sagle, ID
.
Bureau of Reclamation
2013
Downscaled CMIP3 andCMIP5 Climate and Hydrology Projections: Release of Downscaled CMIP5 Projections, Comparison with Preceding Information, and Summary of User Needs
.
U.S. Department of the Interior, Bureau of Reclamation, Technical Services Center
,
Denver, CO
.
Burns
D. A.
Klaus
J.
Hale
M. R.
2007
Recent climate trends and implications for water resources in the catskill mountain region, New York, USA
.
Journal of Hydrology
336
,
155
170
.
https://doi.org/10.1016/j.jhydrol.2006.12.019
.
DRBC
1983
DRBC Resolution 83-13
. .
DRBC
1988
DRBC Resolution 88-22
. .
DRBC
2013
Delaware River Basin Water Code
.
Available from: https://www.state.nj.us/drbc/library/documents/watercode.pdf (accessed 19 February 2019)
.
DRBC
2017
Delaware River Basin Commission Basin Annual Report
.
Available from: https://www.state.nj.us/drbc/about/public/annual-reports.html (accessed 19 February 2019)
.
DRBC
2019a
Delaware River Basin Commission Basin Information
.
Available from: https://www.state.nj.us/drbc/basin/ (accessed 19 February 2019)
.
DRBC
2019b
Dates of DRBC-Declared Drought
. .
Ellis
A. W.
Hawkins
T. W.
Balling
R. C.
Jr.
Gober
P.
2008
Estimating future runoff levels for a semiarid fluvial system in central Arizona, USA
.
Climate Research
35
,
227
239
.
https://doi.org/10.3354/cr00727
.
EPA
2011
National Action Plan: Priorities for Managing Freshwater Resources in a Changing Climate
. .
EPA
2016
National Action Plan Update: Priorities for Managing Freshwater Resources in a Changing Climate
.
Available from: https://acwi.gov/climate_wkg/iwrcc/2016_nap_final_20161129.pdf (accessed 20 February 2019)
.
Frei
A.
Armstrong
R. L.
Clark
M. P.
Serreze
M. C.
2002
Catskill mountain water resources: vulnerability, hydroclimatology, and climate-change sensitivity
.
Annals of the Association of American Geographers
92
,
203
224
.
Friedrichi
K.
Grossman
R. L.
Huntington
J.
Blanken
P. D.
Lenters
J.
Holman
K. D.
Gochis
D.
Livneh
B.
Prairie
J.
Skeie
E.
Healey
N. C.
Dahm
K.
Pearson
C.
Finnessey
T.
Hook
S. J.
Kowalski
T.
2018
Reservoir evaporation in the western United States: current science, challenges, and future needs
.
Bulletin of the American Meteorological Society
99
,
167
187
.
https://doi.org/10.1175/BAMS-D-15-00224.1
.
Haff
L.
Demberger
S.
Baumbach
E.
(eds)
2017
Technical Report for the Delaware Estuary and Basin 2017. Partnership for the Delaware Estuary. PDE Report No. 17-07. 379 pages
.
Hamon
W. R.
1961
Estimating potential evapotranspiration
.
Proceedings of the American Society of Civil Engineers
87
,
107
120
.
Hawkins
T. W.
2015
Simulating the impacts of projected climate change on streamflow hydrology for the chesapeake Bay watershed
.
Annals of the Association of American Geographers
105
(
4
),
527
648
.
https://doi.org/10.1080/00045608.2015.1039108
.
Hawkins
T. W.
Austin
B. J.
2012
Simulating streamflow and the effects of projected climate change on the savage river, Maryland, USA
.
Journal of Water and Climate Change
3
,
28
43
.
https://doi.org/10.2166/wcc.2012.016
.
Intergovernmental Panel on Climate Change (IPCC)
2013
Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change
.
Cambridge University Press
,
Cambridge
,
UK
.
Kang
H.
Sridhar
V.
2017
Assessment of future drought conditions in the chesapeake Bay watershed
.
Journal of the American Water Resources Association
54
(
1
),
160
183
.
https://doi.org/10.1111/1752-1688.12600
.
Mather
J. R.
1978
The Climatic Water Budget in Environmental Analysis
.
Lexington Books
,
Lexington, MA
.
Matonse
A. H.
Pierson
D. C.
Frei
A.
Zion
M. S.
Schneiderman
E. M.
Anandhi
A.
Mukundan
R.
S. M.
2011
Effects of changes in snow pattern and the timing of runoff of NYC water supply system
.
Hydrological Processes
25
,
3278
3288
.
https://doi.org/10.1002/hyp.8121
.
Matonse
A. H.
Pierson
D. C.
Frei
A.
Zion
M. S.
Anandhi
A.
Schneiderman
E.
Wright
B.
2013
Investigating the impact of climate change on New York city's primary water supply
.
Climatic Change
116
,
437
456
.
https://doi.org/10.1007/s10584-012-0515-4
.
Mauer
E. P.
Brekke
L.
Pruitt
T.
Duffy
P. B.
2007
Fine resolution climate projections enhance regional climate change impact studies
.
EOS Transactions of AGU
88
,
504
.
https://doi.org/10.1029/2007EO470006
.
McCabe
G. J.
Ayers
M. A.
1989
Hydrologic effects of climate change in the Delaware River Basin
.
Water Resources Bulletin
25
(
6
),
1231
1242
.
https://doi.org/10.1111/j.1752-1688.1989.tb01335.x
.
McCabe
G. J.
Markstrom
S. L.
2007
A Monthly Water-Balance Model Driven by a Graphical User Interface
.
U.S. Geological Survey Open-File Report 2007–1088
.
U.S. Geological Survey
,
Reston, VA
.
McCabe
G. J.
Wolock
D. M.
1992
Effects of climatic change and climatic variability on the Thornthwaite Moisture Index in the Delaware River Basin
.
Climatic Change
20
,
143
153
.
https://doi.org/10.1007/BF00154172
.
McCabe
G. J.
Wolock
D. M.
1999
Future snowpack conditions in the western United States derived from general circulation model climate simulations
.
Journal of the American Water Resources Association
35
,
1473
1484
.
https://doi.org/10.1111/j.1752-1688.1999.tb04231.x
.
McCabe
G. J.
Wolock
D. M.
2011
Century-scale variability in global annual runoff examined using a water balance model
.
International Journal of Climatology
31
,
1739
1748
.
https://doi.org/10.1002/joc.2198
.
Milly
P. C. D.
Dunne
K. A.
2017
A hydrologic drying bias in water-resource impact analyses of anthropogenic climate change
.
Journal of the American Water Resources Association
53
(
4
),
822
838
.
https://doi.org/10.1111/1752-1688.12538
.
Najjar
R.
1999
The water balance of the Susquehanna River Basin and its response to climate change
.
Journal of Hydrology
219
,
7
19
.
https://doi.org/10.1016/S0022-1694(99)00041-4
.
Najjar
R.
Walker
H. A.
Anderson
P. J.
Barron
E. J.
Bord
R. J.
Gibson
J. R.
Kennedy
V. S.
Knight
C. G.
Megonigal
J. P.
O'Connor
R. E.
Polsky
C. D.
Psuty
N. P.
Richards
B. A.
Sorenson
L. G.
Steele
E. M.
Swanson
R. S.
2000
The potential impacts of climate change on the mid-Atlantic coastal region
.
Climate Research
14
,
219
233
.
Najjar
R.
Patterson
L.
Grahm
S.
2009
Climate simulations of major estuarine watersheds in the mid-Atlantic region of the US
.
Climatic Change
95
,
139
168
.
https://doi.org/10.1007/s10584-008-9521-y
.
Najjar
R. G.
Pyke
C. R.
M. B.
Breitburg
D.
Hershner
C.
Kemp
M.
Howarth
R.
Mulholland
M. R.
Paolisso
M.
Secor
D.
Sellner
K.
Wardrop
D.
Wood
R.
2010
Potential climate-change impacts on the Chesapeake Bay
.
Estuarine, Coastal and Shelf Science
86
,
1
20
.
https://doi.org/10.1016/j.ecss.2009.09.026
.
Nash
J. E.
Sutcliffe
J. V.
1970
River flow forecasting through conceptual models part I: a discussion of principles
.
Journal of Hydrology
10
,
282
290
.
https://doi.org/10.1016/0022-1694(70)90255-6
.
Neff
R.
Chang
H.
Knight
C. G.
Najjar
R. G.
Yarnal
B.
Walker
H. A.
2000
Impact of climate variation and change on mid-Atlantic hydrology and water resources
.
Climate Research
14
,
207
218
.
Ning
L.
Riddle
E. E.
R. S.
2015
Projected changes in climate extremes over the Northeastern United States
.
Journal of Climate
28
,
3289
3310
.
https://doi.org/10.1175/JCLI-D-14-00150.1
.
Priestley
C. H. B.
Taylor
R. J.
1972
On the assessment of surface heat flux and evaporation using large-scale parameters
.
Monthly Weather Review
100
,
81
92
.
https://doi.org/10.1175/1520-0493(1972)100<0081:OTAOSH > 2.3.CO;2
.
Smith
K.
Hawkins
T. W.
2011
Hydroclimate trends in the Susquehanna River Basin
.
The Pennsylvania Geographer
49
(
2
),
34
49
.
Sohl
T. L.
Sayler
K. L.
Bouchard
M. A.
Reker
R. R.
Friesz
A. M.
Bennett
S. L.
Sleeter
B. M.
Sleeter
R. R.
Wilson
T.
Soulard
C.
Knuppe
M.
Van Hofwegen
T.
2014
Spatially explicit modeling of 1992–2100 land cover and forest stand age for the Conterminous United States
.
Ecological Applications
24
(
5
),
1015
1036
.
https://doi.org/10.1890/13-1245.1
.
Sohl
T. L.
Reker
R.
Bouchard
M.
Sayler
K.
Dornbierer
J.
Wika
S.
Quenzer
R.
Friesz
A.
2016
Modeled historical land use and land cover for the conterminous United States
.
Journal of Land Use Science
11
(
4
),
476
499
.
https://doi.org/10.1080/1747423X.2016.1147619
.
Taylor
K. E.
Stouffer
R. J.
Meehl
G. A.
2012
An overview of CMIP5 and the experiment design
.
Bulletin of the American Meteorological Society
93
,
485
498
.
https://doi.org/10.1175/BAMS-D-11-00094.1
.
Thornthwaite
C. W.
1948
An approach towards rational classification of climate
.
The Geographical Review
38
,
55
94
.
USDA
2017
Gridded Soil Survey Geographic (gSSURGO) Database User Guide
. 29
January 29 2017
.
U.S. Department of Agriculture Natural Resources Conservation Services
.
USGS (U.S. Geological Survey)
2015
National Hydrography Dataset (NHDVPlus02)
. .
USGS (U.S. Geological Survey)
2017
Agreement for a Flexible Flow Management Program
.
Available from: https://webapps.usgs.gov/odrm/ffmp/FFMP2017.pdf (accessed 20 February 2019)
.
Williamson
T. N.
Lant
J. G.
Claggett
P. R.
Nystrom
E. A.
Milly
P. C. D.
Nelson
H. L.
Hoffman
S. A.
Colarullo
S. J.
Fischer
J. M.
2015
Summary of Hydrologic Modeling for the Delaware River Basin Using the Water Availability Tool for Environmental Resources (WATER)
.
U.S. Geological Survey Scientific Investigations Report 2015-5143
,
68 p
.
Williamson
T. N.
Nystrom
E. A.
Milly
P. C. D.
2016
Sensitivity of the projected hydroclimatic environment of the Delaware River Basin to formulation of potential evapotranspiration
.
Climatic Change
139
,
215
228
.
https://doi.org/10.1007/s10584-016-1782-2
.
C. J.
2016
Meteorologic, hydrologic, and geographic characterization of historic floods in the Delaware River Basin
.
Middle States Geographer
49
,
63
73
.
C. J.
Hawkins
T. W.
Jantz
C.
Drzyzga
S.
In: Review
impact of changing climate and land cover on flood magnitudes in the Delaware River Basin, USA
.
Journal of the American Water Resources Association
.
(In review.)
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).