Quantifying Hydrological Impacts of Climate Change Uncertainties on a Watershed in Northern Virginia

Forecasted changes to climate were used to model variations in the streamflow characteristics of a northern Virginia catchment. Two emission scenarios were applied from international climate projections along with four general circulation models (GCMs) by using two statistical downscaling methods to drive the hydrological simulations in two future time periods (2046–2065 and 2081–2100). Incorporation of these factors yielded 32 runoff simulation models for a 130-km watershed located in northern Virginia. These models were compared with historical streamflow data from the late 20th century. Changes in streamflow were compared using median, low, and high flows. Results showed a general increase in median flows in both the midand late 21st century. Low flows were projected to decrease, whereas high flows were projected to increase, creating a larger range between low flows and high flows. In addition, statistical tests were conducted to identify the main factors that affected variations in future climate projections. The choice of the downscaling method emerged as the main source of uncertainty. This research quantifies the impacts of climate change as well as uncertainties within climate change projections for regional water resources. Considering the essential role of this watershed for water supply in northern Virginia, the findings of this study illustrate likely impacts of climate change on water supply reliability, supporting climate resiliency in the study area. DOI: 10.1061/(ASCE)HE.19435584.0001860. © 2019 American Society of Civil Engineers. Author keywords: Climate change; Hydrological impact; Scenario development; Water supply; Water resources management.


Introduction
The effects of anthropogenic climate change on the water cycle have been a focus of many studies in recent years (Haddeland et al. 2014;Ichoku and Adegoke 2016;Stocker and Raible 2005). Commonly, recorded time series of observed data and historical events are used to analyze, design, plan, and manage water resource systems for future events. Many of these procedures assume a stationary state for observed data that carry over to future estimates. These assumptions thus do not reflect possible additional needs caused by future climate change and may not be suitable for future plans and designs (Moglen and Rios Vidal 2014). General circulation models (GCMs) are often implemented to predict effects of climate change. GCMs are mathematical representations of climatic variables that incorporate carbon dioxide as one of the major anthropogenic contributors to the greenhouse effect. GCMs use internally coherent physical and, in some cases, chemical characteristics of climate to simulate and quantify the climatic responses to future human-induced activities (Bernstein et al. 2008).
Climate change impact assessments often require a higher resolution than otherwise provided by GCMs. Different downscaling methods have been developed to translate temporal and spatial climate information from coarse-resolution GCM outputs to local and regional scales. Regionalization of the GCMs or downscaling methods can be categorized as either statistical or dynamical procedures (Bernstein et al. 2008). Statistical downscaling methods are based on the premise that the regional climate is conditioned by, first, the large-scale climatic state, and second, by the local physiographic features (Wilby et al. 2004). In this context, quantitative relationships can be established between the large-scale climate variables (predictors) and regional variables (predictands) (Abatzoglou and Brown 2012). Primary types of statistical downscaling methods include (1) weather classification schemes, (2) transfer functions, and (3) weather generators. In contrast, dynamical downscaling methods are based on regional climate models (RCMs) nested within the domain of the GCMs, resulting in higher resolution than the original GCM output. The quality of RCMs is not only dependent upon initial conditions (e.g., soil moisture and soil temperature) but is also dictated by the validity of boundary conditions defined in the GCMs (Wilby et al. 2004). Xu et al. (2005) arrayed a general scheme for assessing the effects of climate change on hydrological regimes based on (1) employment of GCMs representing future climate scenarios as a consequence of increasing greenhouse gases, (2) application of downscaling techniques to downscale GCMs to compatible scales for hydrological modeling, and (3) implementation of hydrological models to simulate the response of hydrological regimes due to climate change. There are several sources of uncertainties in each section of this framework.
As an example, for greenhouse gas (GHG) emission scenarios in the future, the rate of conversion of a specific GHG emission into atmospheric concentrations, the range of responses of various climate models to a given radiative force, and the method of constructing high-resolution information from global climate model outputs all contribute as sources of uncertainty (Mearns et al. 2001). The possible consecutive uncertainties in this framework are known as a "cascade of uncertainties" (Mearns et al. 2001 practice to consider this cascade of uncertainties when developing climate-related scenarios for climate change impact, adaptation, and mitigation assessment studies (Smithson 2002). To incorporate the aforementioned physical uncertainties in the climate system models of this study, combinations of different GHG emission scenarios, GCMs and downscaling methods, and time periods have been considered. Kay et al. (2009) investigated the uncertainty of the climate change impact on flood frequency in England. They considered six different sources of uncertainty: (1) future greenhouse gas emissions, (2) GCMs, (3) downscaling from GCMs (along with a RCM), (4) hydrological models, (5) hydrological model parameters, and (6) internal variability of the climate system. They concluded that the uncertainty arising from future climate models is larger than the corresponding emission scenarios and hydrological modeling, which made GCMs the primary source of uncertainty. Veijalainen et al. (2010) used 20 different GCMs to assess climate change on flooding in Finland. They suggested that GCMs are a greater source of uncertainty than emission scenarios and recommended that outputs of several GCMs be used in climate change impact studies on flooding. Camici et al. (2013) explored the impacts of climate change on flood frequency in central Italy using two GCMs, one emission scenario, and two statistical downscaling methods for the late 21st century. They concluded that the choice of GCMs and downscaling methods play a substantial role in determining the effects of future climate change on flood frequency. Sunyer et al. (2012) used four different regional climate models and five various statistical downscaling methods to estimate rainfall extremes under climate change in Denmark from the period of 2071-2100 using an A1B emission scenario. They observed significant uncertainties in the mean, standard variation, skewness, probability of dry days, and extreme events by using different downscaled projections. They also concluded that a considerable source of uncertainty comes from the choice of climate models. Dobler et al. (2012) quantified various sources of uncertainty in an Alpine watershed using three GCMs, three RCMs, and three bias-correction techniques for downscaling. In their study, uncertainty from bias correction was found to have the highest impact on extreme riverflow projections. They also mentioned that one of the techniques (quantile-quantile mapping) showed more reliable results compared with the other two methods (delta change and local scaling).
In this context, LaFond et al. (2014) downscaled an ensemble of nine GCMs using two different downscaling methods for a watershed in Kentucky and concluded that the delta-change downscaling method was suitable for estimating extreme events in flood impact assessments. Vetter et al. (2017) evaluated the sources of uncertainty of projected climate models for 12 river basins around the globe by applying nine hydrological models, four greenhouse gas concentration trajectories, and five GCMs. They used ANOVA to determine the sources of uncertainties and concluded that the choice of GCMs had the largest share of uncertainty compared with the other factors.
In recent years, different aspects of streamflow characteristics have been studied, with the focus on the impacts of climate change. Percent change in the magnitude of runoff is the most common streamflow attribute that is reported. For instance, Acharya (2017) investigated the impact of climate change on streamflow in northern Alabama. That study compared average change in projected future monthly and total annual streamflow with historical observations. Although monthly changes varied from þ40% to −48% under different scenarios, a general decreasing streamflow pattern in the watershed was observed.
Another common reported change in streamflow characteristics is the timing. For example, Eisner et al. (2017) explored seasonal streamflow across 11 river basins across the globe under climate change. Overall, they found only a minor impact on changes in the timing of high-flow and low-flow periods. However, climate change was shown to have a considerable effect on the magnitude of seasonal streamflow in two basins.
In this study, the authors explore temporal measures of streamflow magnitude as well as changes in the frequency of high and low flows. The study watershed is located in the Virginia suburbs of the Washington, DC, metropolitan area. Recent observations and previous studies in this region (US mid-Atlantic) have indicated a change in precipitation trends and an increase in air temperature (Bernstein et al. 2008;Moglen and Rios Vidal 2014). Maldonado and Moglen (2012) studied the effects of climate change and landuse change on the Occoquan Reservoir with a focus on low-flow variations. They concluded that land-use/demand change has a more substantial impact than climate change on the water supply system, particularly in the mid-21st century of the study area. Moglen and Rios Vidal (2014) investigated the effects of climate change on stormwater infrastructure in the mid-Atlantic region of the United States. They indicated that detention ponds, and likely other stormwater infrastructure systems, would not function as previously envisioned due to the changes in future precipitation intensity and volume. Observed trends of climate change in the study area, as well as this watershed's important role in providing water supply in northern Virginia, necessitates a more comprehensive climate change impact study.
To account for different sources of uncertainties, this work employs future precipitation and surface air temperature time series from four GCMs using two GHG emission scenarios. For regional use, these time series are subsequently downscaled to the watershed by applying two statistical downscaling methods for the projection periods of 2046-2065 and 2081-2100. The objectives of this study are to (1) simulate future streamflow conditions in the study area for the mid-and late 21st century, (2) analyze climate change impacts on the watershed quantifying flow regime alterations, and (3) explore the main sources of uncertainty that arise from differing emission scenarios, GCMs, downscaling methods, and projection periods.

Study Area and Data
The study area is a part of the mid-Atlantic region of the United States and has four distinctive seasons. The Broad Run watershed with a drainage area of 130 km 2 is part of the larger Occoquan watershed in northern Virginia, located immediately west of Washington, DC (Fig. 1). Broad Run is a major contributor to the Lake Manassas water supply. Lake Manassas is a primary drinkingwater source for the City of Manassas, City of Manassas Park, and the surrounding areas (Liu 2011).
Rapid population growth in northern Virginia has increased the demand for water supply, calling for new strategies to meet future needs. According to the US Census Bureau (2019), population in one of the major counties within the Occoquan watershed (Prince William County) increased 260% from 1970 to 2010. Urbanization is an ongoing theme in the study area, which can greatly influence a watershed's response to precipitation (Leopold 1968). The Broad Run watershed has moderate to high slopes with elevations ranging from 122 to 305 m. The general structure of soils is loamy, with approximately 67% in Group B of the Natural Resources Conservation Service (NRCS) hydrological soil groups (Liu 2011).
Daily meteorological data are collected at several weather stations within and near the study location. Cloud coverage, dew-point temperatures, solar radiation, and wind speeds are obtained from Washington Dulles International Airport. Air temperature and precipitation data are obtained from three stations within the Broad Run watershed. Observed annual precipitation varies from 670 to 1,360 mm with a yearly average of 960 mm. Observed mean annual surface air temperature is roughly 12.4°C with a maximum monthly mean surface air temperature of 30.6°C in July and minimum monthly mean surface air temperature of −4.5°C in January. Table 1 depicts the observed meteorological data for the study area.
Evapotranspiration was simulated using Hamon (1963)'s potential evapotranspiration (PET) formula, which relates potential evapotranspiration to maximum possible incoming radiant energy and saturation vapor pressure (Lu et al. 2005). The potential moisture-holding capacity of the air was calculated at the prevailing air temperature. Radiant energy was computed using daylight hours, sunset hour angle, Julian day of the year, and latitude of the study area.

Modeling Climate Change Impact
Climate change impact studies often consider an ensemble of models to address the uncertainties arising from each factor of the modeling process (Camici et al. 2013;Dobler et al. 2012;Moglen and Rios Vidal 2014). A similar procedure has been implemented  in this study. For this purpose, two emission scenarios, four GCMs, two downscaling methods, and two simulation periods were used to develop future time series for precipitation (P) and surface air temperatures (T). In addition, PET time series were calculated using the method of Hamon (1963), in which projected maximum and minimum surface air temperatures were used. The process of developing projected maximum and minimum surface air temperatures were analogous to the method used to develop surface air temperature time series. Cloud coverage, dew-point temperatures, solar radiation, wind speeds, and land use were assumed to be stationary for the two future time periods (2046-2065 and 2081-2100) investigated. Possible changes in soil type and infiltration rate, stream morphology, deep groundwater recharge, and losses are assumed to remain the same in future projections.

Emission Scenarios
GHG emissions are exacerbated by anthropogenic activity. Quantifying the evolution of future trends of GHG emissions is complicated. The driving forces of GHG emission scenarios are highly dynamic, connecting several components such as demographic development, socioeconomic progress, agricultural practices, and technological changes (Smithson 2002). Using available data and existing knowledge of GHG patterns, the Intergovernmental Panel on Climate Change (IPCC) has developed different storylines and scenario families to cover the extensive range of possible uncertainties and future changes (Bernstein et al. 2008). Four scenario families have been established as the main driving forces of GHG emissions (Bernstein et al. 2008). To address the uncertainty within GHG emission scenarios in this study, two extremes of the Special Report on Emissions Scenarios (SRES) (IPCC 2007) spectrum have been considered, giving rise to the full range of emission possibilities. On the upper side of the range, the A2 scenario focuses on regional and economic growth, accounting for increasing global population, locally oriented economic growth, and fragmented and slower technological advancement. On the lower side of the range, the B1 scenario assumes greater environmental awareness and global growth attributed to convergent global population after the mid-21st century, introduction of clean technologies, and environmental sustainability (Bernstein et al. 2008).

GCMs
Each GCM represents the Earth's climate system under different initial and boundary conditions. To address the uncertainty specifically arising from the GCMs, an ensemble of four GCMs have been considered. However, only one realization of each GCM was employed, meaning that the internal variability simulated in each GCM is not reflected ( Table 2).

Statistical Downscaling of Precipitation and Surface Air Temperature Time Series
Simulated data from GCMs are primarily used to demonstrate the Earth's climate system on a global scale (e.g., 1.875°× 1.875°). However, for regional use, these models should be translated to the scale of the watershed being studied. For this reason, the use of an appropriate downscaling method is required. In this study, two statistical downscaling methods, delta change (Anandhi et al. 2011;Arnell 1996;Camici et al. 2013;Wilby et al. 2004) and quantile mapping (Bennett et al. 2014;Camici et al. 2013;Maurer and Hidalgo 2008;Piani et al. 2010;Thrasher et al. 2012;Wood et al. 2004), were applied to produce the needed regional data sets. The essential advantages of using statistical downscaling are computational efficiency and the ability to develop GCM ensembles from multiple GCMs (Abatzoglou and Brown 2012). Moreover, this method can directly assimilate site-specific observations, which can be crucial for many operative and managerial climate change impacts studies (Wilby et al. 2004).

Delta Change Method
The delta change (DC) method is one of the simplest forms of downscaling transfer functions and change factor methodologies.
In this method, different formulations based on temporal scales and mathematical operations (generally additive or multiplicative) are implemented. Many researchers use additive formulations as an estimate of the absolute change in temperature and the relative change for future precipitation time series (Anandhi et al. 2011).
In the case of downscaling, surface air temperature data using additive change factors may be preferable because multiplicative change factors may encounter an undefined or unrealistically small or large variable. On the other hand, some meteorological variables, such as precipitation, have a lower limit (i.e., 0), and additive change factors must be checked for negative values. In this study, an additive change factor [Eqs.
(1) and (2)] was applied to produce future surface air temperatures projections and a multiplicative change factor [Eqs. (3) and (4)] was applied to create future precipitation time series where FC at = future climate surface air temperature; HC at = historical climate surface air temperature; δ at = change factor for surface air temperatures; LH at = local historical surface air temperatures; and LF at = local future surface air temperature where FC pr = future climate precipitation; HC pr = historical climate precipitation; δ pr = change factor for precipitation; LH pr = local historical precipitation; and LF pr = local future precipitation.

Quantile Mapping Method
In the quantile mapping (QM) method, cumulative distribution functions (CDFs) of observed and GCMs time series are generated for the variables of interest (i.e., temperature and precipitation). One CDF is based on the GCMs simulations and the other is based on the aggregated observations. The mapping procedure is a simple nonparametric lookup that matches statistical moments of the observed data to the simulated ones. The QM method can reflect changes in the mean and variance of climate variables matching all statistical moments of GCMs with the corresponding observed variables, whereas the delta change method can only consider the mean change of the GCMs variables. Different variations of QM methods have been developed by researchers.
In this study, bijective (predictors are the same parameters as the predictands) and parameter-free (empirical cumulative distribution function) QM were implemented. The equidistant cumulative distribution function (EDCDF) matching method was used for downscaling surface air temperatures [Eq. (5)] (Li et al. 2010). EDCDF has limitations for bias-correcting precipitation where negative values are not acceptable. As a result, some researchers have tried to fit a theoretical distribution to their data to overcome this shortcoming. Alternatively, in this study, the equiratio cumulative distribution function (ERCDF) matching method was implemented to downscale precipitation data for prevailing this drawback [Eq. (6)] (Wang and Chen 2014) Simulation Time Period In this study, impacts of future climate change are analyzed over two time periods, the mid-21st century and late 21st century, representing the near and far future, respectively. The period of the mid-21st century consists of simulations from 2046 to 2065, and the period of the late 21st century entails simulations from 2081 to 2100. Table 3 summarizes the historical and GCM data examined in this study.

Hydrological Simulation Model
Reliable rainfall-runoff modeling for a watershed is a priority for measuring streamflow alteration in climate change impact studies. In this paper, hydrological simulations were performed using the Hydrological Simulation Program-FORTRAN version 12.2 (HSPF). HSPF is jointly supported by the United States Environmental Protection Agency and US Geological Survey. HSPF is a comprehensive process-based watershed model that has the ability to integrate watershed hydrology and water quality simulations. The extended period structure of HSPF allows users to simulate pervious and impervious land surfaces as well as in-stream hydraulic and sediment-chemical interactions (USEPA 2019). HSPF has been widely used as an analytical tool for planning, designing, and operating water resources systems throughout North America and numerous countries around the globe with different climatic regimes (USEPA 2019).
For this watershed, delineation and segmentation have been performed based on local geographic characteristics such as topography, land use, and soil properties (Liu 2011). The HSPF model contains three application modules: pervious land segments (PERLND), impervious land segments (IMPLND), and reaches (RCHRES). These modules can be used to simulate runoff and water quality constituents in pervious and impervious lands as well as free-flow reaches/well-mixed impoundments. In this study, the HSPF model version 12.2 developed by Liu (2011) was used. This model was calibrated for a 4-year period (2002)(2003)(2004)(2005) and validated for a 2-year period (2006)(2007). Liu (2011) used several parameters to calibrate the hydrological model using observed data from stream station ST70 (Fig. 1). This model was calibrated by maximizing the Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe 1970) and by minimizing the percent bias (PBIAS) (Moriasi et al. 2007). The calibrated model has been found be quite robust and able to successfully simulate a range of scenarios for this watershed including flows during wet periods.

Flow-Regime Indicators
Despite the development of several streamflow alteration indicators in hydrological analyses in previous studies, in many cases the impacts of future climate change on runoff have been expressed solely by using percent change of the magnitude of runoff with respect to natural conditions in the baseline period. Olden and Poff (2003) reviewed 171 different hydrological indexes for characterizing streamflow regimes. The main classifications of indexes in their study were magnitude, frequency, duration, timing, and rate of change in flow events. These indexes are calculated daily, monthly, or annually. Researchers have used and adapted various indexes to address the alterations in streamflow due to climate change.
Döll and Schmied (2012) studied the impact of climate change on freshwater resources using mean annual runoff (MAR), statistical low and high flows (monthly discharges that exceed the 90th and 10th percentile of flows), and mean seasonal discharge. Laizé et al.
(2014) adopted a set of monthly flow-regimes indicators (MFRIs) based on indicators of hydrologic alteration (IHA) (Richter et al. 1996) for assessing flow alterations and ecological risks for pan-European rivers. In this context, Thompson et al. (2014) used MFRIs to assess environmental flows and ecological risks for the Mekong River in southeast Asia. In the present study, a set of indexes to quantify streamflow alterations was adopted as provided in Table 4.

Surface Air Temperature
All models predicted an annual increase in average surface air temperature ranging from 0.9°C to 4.8°C. The smallest change in annual surface air temperature was from the MRI model using Emission scenario B1 and the DC downscaling method for the mid-21st century horizon. The highest annual surface air temperature was projected by the GFDL model using Emission scenario A2 and the QM downscaling method in the late 21st century. On average, an increase of 11% and 19% for annual surface air temperature was predicted across all models for the 2046-2065 and 2081-2100 periods, respectively. Annual maximum surface air temperature was also projected to increase by 7% and 13% on average for the first and second time periods. The highest change in annual maximum surface air temperature was 24%, which was projected by the GFDL model using Emission scenario A2 and the QM downscaling method for the second time period. The smallest change in annual maximum surface air temperature was 4%, projected by the CSIRO model using Emission scenario B1 and the DC downscaling method for the mid-21st century. Annual minimum surface air temperatures showed a 21% and 36% increase for the first and second time periods, respectively, compared to the annual minimum surface air temperatures in the baseline period.
All models showed increases in maximum and minimum surface air temperatures, with minimum temperatures predicted to rise at slightly higher amounts. This led to a decrease in the annual average of diurnal temperature range, which is the difference between the daily maximum and minimum temperatures. This range decreased more in the second period. Increases in maximum and minimum temperatures could affect future physiological processes of plant ecosystems and biogeochemical cycling modes in significant ways, including phenophase transitions to spring and budburst timing, autumn leaf abscission, and cold hardening (BACC Author Team 2008). Fig. 2 illustrates monthly surface air temperatures for all models in both the mid-and late 21st century periods. The highest rise (9%) in average monthly surface air temperature was projected during the summer months (June, July, and August). August had peak monthly average surface air temperatures of 27.3°C and 31.6°C for the mid-and late 21st century, respectively, using the GFDL model, Emission scenario A2, and the QM downscaling method. Maximum monthly surface air temperature, like the average monthly surface air temperatures, increase most during the summer months. However, monthly average minimum surface air temperature was anticipated to change more in the autumn (September, October, and November) with an increase of 1.7°C and 2.8°C for the first and second future periods, respectively.

Precipitation
Mean annual precipitation was projected to increase, ranging from 2% to 17%, with the medians of 9% and 11% increases for the midand late 21st century, respectively. Among the GCMs, CSIRO showed the greatest change for the first time period, whereas GFDL showed the highest change for the second period. For both periods, MRI projected the least amount of increase in precipitation. Overall, the A2 emission scenario forecasted more precipitation than the B1 scenario. The DC downscaling method produced slightly higher increments of precipitation than the QM method. The median precipitation volumes were higher for all months, except September, with the highest projected precipitation for May in both periods. The seasonal increases in median precipitation ranged from 3% in autumn to 13% in spring. Fig. 3 displays the monthly variations of precipitation for the mid-and late 21st centuries.   monotonically increasing density grades were applied to show the composition and proportion of summed densities for each period. Increases in the median annual precipitation and the median annual mean surface air temperature were projected for both periods. The median of annual mean surface air temperatures was expected to increase 1.6°C and 2.3°C, for mid-and late 21st century periods, respectively. Changes in the median annual precipitation were projected to increase 9% and 11% for the mid-and late 21st century periods, respectively. Overall, the ranges of both median annual precipitation and median annual surface air temperature of projected models were higher for the late 21st century period.

Streamflow Alterations
A statistical test was applied to explore the differences between the projected mean annual flows and the historical data. The null hypothesis was that the mean annual flows in the future would be the same as those experienced in the late 20th century. In more detail, the null hypothesis assumed that the mean annual flows in the projected models came from the same population as the observed data. A nonparametric method, the Wilcoxon signed rank test (Wilcoxon et al. 1963), was implemented. Projection models in which the null hypothesis with 95% confidence interval (α ¼ 0.05) was rejected are presented in Table 5.
Rejecting the null hypotheses for these cases indicated statistically significant changes in mean annual flows with respect to historical values. The Hodges-Lehmann estimator (Hodges and Lehmann 1983) was subsequently implemented to estimate the magnitude of differences between the projected models and the historical data. As presented in Table 5, some of the projections showed up to a 30% change in the magnitude of mean annual flows.   Table 6 provides the details of changes in the characteristics of flow in the study area. In general, model projections showed increasing trends in the mean annual flows. These results are in agreement with previous findings, which signified an increase in mean monthly and yearly flows (Stagge and Moglen 2017). Projections using the DC method showed higher changes in the mid-21st century, whereas the QM method projected greater changes in late 21st century. Moreover, changes using A2 scenarios were more pronounced than B1 scenarios.
Q 10 is an indicator of high flows, which is defined as the flow that exceeds 90% of the flows on a monthly basis in the overall period. Percent change in Q 10 for any given model is calculated as the percent change between the Q 10 of the projected model and the historical data. The same procedure was performed for Q 90 and Q IQR , as indicated in Table 6. Unlike Model group QM-B1, all other model groups showed increases in high flows (Q 10 ) for the future periods. On the other hand, low flows showed decreasing trends, suggesting that wider intervals between higher and lower deciles result in higher variabilities in the future flows. In other words, the model projected a higher amount of high flows and a higher amount of low flows in the future. An increase in the amount of high flows will amplify the risk of floods, and a decrease in the amount of low flows will intensify the dry spells and cause ecological risks. As presented in Table 6, this trend also affected the interquartile range. The interquartile range indicates the dispersion of the central 50% of data. On average, this range increased by 13% and 14% for the mid-and late 21st century periods, respectively.
Percent change in the number of months above and below the threshold was calculated to determine not just the change in magnitude but also change in frequency of the months that have higher or lower flows. To compute percent change in number of months above and below the threshold, the numbers of months that exceeded the Q 10 and Q 90 thresholds in the historical data were first calculated. Then, the numbers of monthly flows in the projection models that exceeded historical Q 10 and Q 90 thresholds were counted. The percent change between these numbers is given in Table 6 for each of the different group models during the midand late 21st centuries. As presented in Table 6, generally, the number of months that exceed the Q 10 and Q 90 thresholds increase in both periods. As a result, climate change can cause a shift in location, an increase in scale, and also an increase in frequency of high and low flows in this catchment.
To quantify the relative impacts on the different projection models examined, an ANOVA model was employed to determine which modeling factor had the largest impact on model outcomes. Two different data transformation methods were used prior to the ANOVA design. The first data transformation was applied using a power transformation (Helsel and Hirsch 2002) to provide an approximately normal and constant variance data set. Because water flow regimes are usually positively skewed, the log transformation was applied as the first data transformation method. The rank transformation method (Helsel and Hirsch 2002) was performed as the second data transformation method with the benefit of being an asymptotically distribution-free method. The nonparametric fourway ANOVA with fixed factors was thus designed in consideration of four different impacting factors [i.e., choice of emission scenario (two levels), choice of GCM (four levels), choice of downscaling method (two levels), and choice of simulation period (two levels)]. In this design, there were four main effects and 11 interactions among factors.
In practice, ANOVA methods are intended to help parse experimental data. This case is an example of using ANOVA for simulation data identifying the degree of variability in the models. In principle, ANOVA was used to determine the relative effects of different factors in future climate projections. In this design, Type III sum of squares was used (Langsrud 2003). However, because the design was orthogonal, meaning that there were equal amounts of simulations in every model, the type of sum of squares did not affect the overall results. The null hypothesis for the rank transformation method, as one of the forms of data transformation, looked for heterogeneity of the mean rank as an estimate of the median among all model groups. The null hypothesis was determined as if the median was the same between the model groups for all factor levels.
The results of ANOVA using both data transformation methods presented similar outcomes, both indicating a rejection of the null hypothesis. The rejection of the null hypothesis gives the conclusion that the medians of the mean annual flows differed between model groups. Among the main factors, the downscaling method with α-value denoted that responses for levels DC or QM were significantly different. Effects of other levels for the remainder of the main factors were not statistically significant, suggesting that only the choice of the downscaling method, among other sets of choices, had a statistically significant effect on changes in the responses (mean annual flows).  Fig. 5 shows that there were some variations across the interactions of the factors. However, these variations were not statistically significant, indicating no evidence of two-way or more than twoway interactions between factor levels.

Conclusions
The hydrological impacts of climate change on a water supply watershed located in northern Virginia were examined in this paper. An ensemble of model groups was developed consisting of two emission scenarios, four GCMs, two downscaling methods, and two periods, resulting in a total of 32 model projections. There were two classes of findings in this study: changes in physical behavior (temperature, precipitation, and streamflow) and an investigation of the sensitivity of the resulting streamflow projections on the parameters of the model grouping used.
With regards to physical behavior, all projection models showed increases in mean annual surface air temperatures. The shift in location and dispersal in mean annual surface air temperatures were higher in the late 21st century than the mid-21st century. Also, both the maximum and minimum surface air temperatures were projected to increase in both time periods. Overall, these changes are projected to affect the lower tropospheric temperature profile because higher amounts of precipitation are anticipated in the study area in both future periods. The precipitation depth was predicted to be slightly higher in the late 21st century by about 2%. The higher temperatures in this period resulted in a 5% increase in evapotranspiration and a 3% reduction in water yield compared with the mid-21st century.
Changes in the aforementioned meteorological variables led to changes in streamflow behavior. To compute the alteration in streamflow regimes, a nonparametric test was carried out to compare differences in mean annual flows. Some of the model projections estimated up to a 30% positive shift in the median of mean annual flows. On average, an 8% change in the median of mean annual flows was anticipated for both the mid-and late 21st century. It is expected that high flows will increase in both the mid-and late 21st century, and the decrease of low flows in the late 21st century will be greater than during the midcentury. At the same time, Q IQR , which is understood to represent the dispersion of the scale parameter, was projected to increase in both time horizons. These changes were more prominent using the QM downscaling method. In essence, these dynamics can result in projected flows with higher intensity runoff, resulting changes in the statistical moments of the flood probability distribution function.
Using an ANOVA test, the sensitivity of mean annual streamflow to the method used for downscaling was found to be statistically significant. Other sensitivities of the modeling approach (emission scenario, GCM choice, and simulation period) were investigated but were not found to result in a statistical difference in simulated mean annual streamflow.
There are several important implications that may be drawn from this study. First, climate forecasts provide a useful tool for estimating future water supply conditions. Forecasted changes in temperature and precipitation, not surprisingly, result in changes in streamflow and thus, water supply as described extensively in the preceding sections. Second, this study shows mixed implications for actual water supply in the Occoquan system with projected increases in mean annual streamflow (i.e., stronger supply) but accompanied by greater variability (i.e., larger uncertainty) in that supply. Finally, there are implications for the research community in that different downscaling methods produce statistically different streamflow forecasts. This indicates a need for more study to determine the most appropriate downscaling approaches for support of water supply planning and decision making.