4>EPA
United States
Environmental Protection
Agency
Final Report
Impact of Urbanization on the
Hydrology of the Pocono Creek
Watershed: A Model Study
-------
EPA/600/R-07/006
September 2006
Impact of Urbanization on the Hydrology
of Pocono Creek Watershed: A Model
Study
By
Mohamed M. Hantush
U.S. Environmental Protection Agency
Cincinnati, Ohio 45268
Latif Kalin
Oak Ridge Institute for Science and Education (ORISE)
U.S. Environmental Protection Agency
Cincinnati, Ohio 45268
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, Ohio 45268
-------
Notice
The U.S. Environmental Protection Agency through its Office of Research and
Development funded the research described here. It has been subjected to the Agency's peer
and administrative review and has been approved for publication as an EPA document.
This research was supported in part by an appointment to the Post Doctoral Research
Program at the National Risk Management Research Laboratory, administered by the Oak
Ridge Institute for Science and Education through Interagency Agreement No. DW89939836
between the U.S. Department of Energy and the U.S. Environmental Protection Agency.
Mention of trade names or commercial products does not constitute endorsement or
recommendation for use.
All research projects making conclusions or recommendations based on environmentally
related measurements and funded by the Environmental Protection Agency are required to
comply with the Agency Quality Assurance Program. This report has been subjected to
QA/QC review. The report presented a mathematical framework for modeling hydrology of a
watershed and did not involve collection and analysis of environmental measurements.
-------
Foreword
The U.S. Environmental Protection Agency (EPA) is charged by Congress with protecting
the Nation's land, air, and water resources. Under a mandate of national environmental laws,
the Agency strives to formulate and implement actions leading to a compatible balance
between human activities and the ability of natural systems to support and nurture life. To
meet this mandate, EPA's research program is providing data and technical support for
solving environmental problems today and building a science knowledge base necessary to
manage our ecological resources wisely, understand how pollutants affect our health, and
prevent or reduce environmental risks in the future.
The National Risk Management Research Laboratory (NRMRL) is the Agency's center for
investigation of technological and management approaches for preventing and reducing risks
from pollution that threaten human health and the environment. The focus of the Laboratory's
research program is on methods and their cost-effectiveness for prevention and control of
pollution to air, land, water, and subsurface resources; protection of water quality in public
water systems; remediation of contaminated sites, sediments and ground water; prevention
and control of indoor air pollution; and restoration of ecosystems. NRMRL collaborates with
both public and private sector partners to foster technologies that reduce the cost of
compliance and to anticipate emerging problems. NRMRL's research provides solutions to
environmental problems by: developing and promoting technologies that protect and improve
the environment; advancing scientific and engineering information to support regulatory and
policy decisions; and providing the technical support and information transfer to ensure
implementation of environmental regulations and strategies at the national, state, and
community levels.
This publication has been produced as part of the Laboratory's strategic long-term research
plan. It is published and made available by EPA's Office of Research and Development to
assist the user community and to link researchers with their clients.
Sally Gutierrez, Director
National Risk Management Research Laboratory
11
-------
Abstract
The Pocono Creek watershed located in Monroe County, PA, is threatened by high
population growth and urbanization. Of concern specifically is the potential impact of future
developments in the watershed on the reduction of base flow and the consequent risk of
degradation of wild brown trout habitats in Pocono Creek. Anticipated increase in
imperviousness, on the other hand, is expected to elevate flood risk and the associated
environmental damage. A watershed hydrology based modeling study was initiated by the
U.S. EPA in collaboration with the U.S. Geological Survey and the Pennsylvania Fish and
Boat Commission to assist Monroe County in planning for sustainable future development in
the Pocono Creek watershed.
The Soil and Water Assessment Tool (SWAT) is selected to model the impact of
projected future build out in the Pocono Creek watershed on the hydrologic response thereof.
The model is successfully calibrated and validated for two sources of precipitation data,
raingauge and Next Generation Weather Radar_(NEXRAD) hourly precipitation data. The
results clearly show that NEXRAD is an effective and economic alternative source of spatio-
temporal precipitation, and that future modeling studies in ungauged watersheds may benefit
from the use of NEXRAD rainfall data.
Ensemble model forecast is constructed using time series analysis and Monte Carlo (MC)
simulations to evaluate model predictive uncertainty. The MC simulations over a 20-year
long period yielded an ensemble of rating curves of which the median and 95% confidence
band of daily streamflows are plotted. These plots allow for the construction of the 95%
confidence band for design flows corresponding to any given recurrence or return period.
SWAT simulated daily streamflow rates in the range 2 to 11 (mVs) show the least
uncertainty. Computed daily streamflow rates below 2 m3/s have the greatest uncertainty,
whereas for flows higher than 11 m3/s uncertainty is moderate.
MC simulations over a 20-year period predict that, on the average, daily base flow would
be reduced by 31% based on the projected build out in the watershed. However, predicted
change in the average of daily streamflow (averaged over the entire simulation period)
appears to be negligible. The computed low-flow index, 7Q10, is expected to decline by
11%, and the monthly median daily flow is expected to be reduced by 10% on the average.
The monthly peak of simulated daily flows and annual maximum daily flow on the average
are predicted to increase by 21% and 19%, respectively. Watershed-averaged groundwater
recharge is predicted to decline by 31% due to the projected land use changes. The median of
the MC simulated flow duration curves shows that in general the likelihood that the
watershed will experience high and low streamflows will increase with the projected
urbanization.
An index methodology is developed to rank seven subwatersheds composing the modeled
portion of the Pocono Creek watershed based on their relative impact on watershed response
to anticipated land developments. The first index, a, signifies the absolute impact of a
particular catchment area on the watershed response. The second index, /?, is a normalized
in
-------
by the percentage area of the sub-catchment, and therefore describes the impact per area of
land use changes. With a few exceptions, a and |3 indices produce similar rankings among
the 7 catchment areas for 7Q10, monthly median of daily flow, and annual maximum daily
flow. These ranking results may be related to groundwater recharge, area, topographic
features, and proximity to the streamflow gauge station. The very downstream catchment
area 7 ranked first in terms of impact on annual maximum daily flows, and second in terms
of impact on 7Q10 and monthly median daily flows. Catchment area 4 associated with the
highest groundwater recharge was ranked first and second for impact on 7Q10 based on a
and |3 indices, respectively. Areas characterized by steep topography and intense wetlands
ranked low, some times the lowest, with respect to impact on the three design flows.
The results of this model study point toward significant changes in low as well as high
flow regimes, should the Pocono Creek watershed experience land use changes consistent
with the projected build out in the watershed. Management measures may be taken in the
future to minimize the predicted changes in the watershed hydrology.
Keywords: watershed, radar, rain gauges, hydrologic modeling, distributed parameter,
SWAT, critical source area, land use change, uncertainty, forecast, Monte Carlo,
time series analysis, Pocono Creek, urbanization, base flow
IV
-------
Table of Contents
Notice i
Foreword ii
Abstract iii
Figures vii
Tables ix
Acknowledgments x
1 Introduction 1
1.1 Pocono Creek Watershed 1
1.2 Objectives 2
1.3 Report Organization 3
1.4 Summary 4
2 Model Selection 5
2.1 Background 5
2.2 Summary of Hydrologic Component of SWAT Model 6
2.3 Summary 9
3 Watershed Model 10
3.1 Study Area 10
3.2 Base flow Separation 11
3.3 Next Generation Weather Radar (NEXRAD) 12
3.4 Model Calibration 12
3.5 Model Validation 18
3.6 Conclusions 20
4 Model Prediction Uncertainty 22
4.1 Introduction 22
4.2 Structural and Model Parameters Uncertainty 23
4.2.1 Prediction Error 23
4.2.2 Nonparametric Probability Distribution 28
4.2.3 Nonparametric Random Generation 29
4.3 Model Forecast Evaluation 30
4.4 Uncertainty Due to Precipitation Variability 33
4.5 Monte Carlo Simulations 34
4.6 Conclusions 37
5 Impact of Land Use Changes 39
5.1 Introduction 39
5.2 Present Land Use Map and Future Build Out 39
5.3 Monte Carlo Simulation 41
5.4 Predicted Changes in Streamflow 41
5.5 Predicted Changes to Flow Frequency and Duration 45
5.6 Conclusions 45
v
-------
6 Critical Source Areas 47
6.1 Introduction 47
6.2 Methodology 47
6.3 Results 48
6.4 Conclusions 51
7 Summary and Conclusions 52
8 References 54
VI
-------
Figures
Figure 1. Pocono Creek Watershed: Location, topography, USGS stream gauge,
climate stations and NEXRAD cell centroids (adapted from Kalin and
Hantush, 2006a) 2
Figure 2. Pocono Creek watershed: a) Land use map; b) Soil map (STATSGO) 12
Figure 3. Calibration procedure in the SWAT model (adapted from Santhi et al.,
2001) 14
Figure 4. Comparison of NEXRAD estimated precipitation with measurements
from local gauge stations. Values shown are in mm 15
Figure 5. Measured stream flow (SF) and estimated surface runoff (SR) plotted
against SWAT simulated counterparts for the calibration (7/1/2002-
5/31/2004) period. NEXRAD estimates are used as precipitation data
source. Top two panels show monthly results, and the bottom panel
shows daily results 17
Figure 6. Measured streamflow (SF) and estimated surface runoff (SR) plotted
against SWAT simulated counterparts for the validation period
(6/1/2004-4/30/2005). NEXRAD estimates are used as precipitation data
source. Top and bottom panels denote monthly and daily results,
respectively 18
Figure 7. Measured daily streamflows plotted against SWAT simulated
counterparts for the post-validation period (5/1/2005-9/30/2005).
NEXRAD estimates are used as precipitation data source 20
Figure 8. Time series of model prediction errors (Observed - SWAT simulated) in
the calibration and validation time period (7/1/2002-4/30/2005): (a)
daily simulations, (b) 3-day moving average 24
Figure 9. Autocorrelation (ACF) and partial autocorrelation (PACF) functions of:
(a) prediction error sh (b) first difference of residual error Vst = st - st.i,
(c) residuals (wr) for ARIMA(1,1,1) x (0,0,1)3. The lag is in units of
days 26
Figure 10. Frequency histogram, Gaussian PDF, and the nonparametric model fit to
the computed series of residuals w^. Both the histogram of residuals and
normal distribution has zero mean and standard deviationcr^ 29
Figure 11. SWAT model ensemble forecast (median and 95% confidence band) and
measured streamflows during the validation period (12/16/2004-
4/30/2005). A total of 9 measurements lie outside the 95% confidence
band (6.7 %). The simulated and measured values are 3-day averages 32
Figure 12. SWAT model forecast (median and 95% confidence band) and measured
streamflows during the post-validation period (5/1/2005-9/30/2005). The
simulated and measured values are 3-day averages 32
Figure 13. SWAT model 3-day average ensemble forecast (median and 95%
confidence band) and measured daily streamflows during the validation
period (12/16/2004-4/30/2005). A total of 15 measurements fall outside
the 95% confidence band (11 %) 33
vn
-------
Figure 14. Median and 95% confidence band for the MC simulated duration curve.
The flow duration curves are generated from 500 replicas of 20 years
daily streamflows, ot 37
Figure 15. Coefficient of variation computed from ensemble of flow duration
curves versus probability of exceedance in the left panel P(X > x), and
cumulative probability P(X
-------
Tables
Table 1. Data used in the model construct and their sources 11
Table 2. Mass Balance Error (MBE), coefficient of determination, R2, and Nash-
Sutcliffe efficiency, ENS for Stream Flow (SF), Base flow (BF) and Surface
Runoff (SR) during the calibration period with NEXRAD as precipitation
data source. Values in parenthesis indicate values when rain gauges are used
as precipitation data source 17
Table 3. Mass Balance Error (MBE), coefficient of determination, R2, and Nash-
Sutcliffe efficiency, ENS for Stream Flow (SF), Base flow (BF) and Surface
Runoff (SR) during the validation period, g and n represent gauge and
NEXRAD input precipitations, respectively 19
Table 4. Values of residual variance, Final prediction error (FPF), Akaike's
information criterion (AIC), and Bayesian information criterion (BIC) for
various models applied to model prediction error data, et 27
Table 5. Summary statistics of ensemble streamflow forecast based on most recent
land use (Figure 2a). Values in the table are obtained from 500 time series of
synthesized daily streamflows 35
Table 6. Summary statistics of computed streamflow characteristics at the USGS
gauge station. The results are derived from 50 MC simulations each 20
years long 42
Table 7. Computed indices a and /? for the 7 management areas for 7Q10, monthly
median of daily flows and annual maximum daily flows 49
IX
-------
Acknowledgments
The U.S. Environmental Protection Agency, through its Office of Research and
Development, funded the research described here through in-house efforts and in part by an
appointment to the Postgraduate Research Program at the National Risk Management
Research Laboratory administered by the Oak Ridge Institute for Science and Education
through an interagency agreement between the U.S. Department of Energy and the U.S.
Environmental Protection Agency. The report benefited from the constructive review
comments of Dennis Lai and Rao S. Govindaraju.
x
-------
1 Introduction
1.1 Pocono Creek Watershed
Pocono Creek watershed is a 46.5 square miles (120 km2) basin located in Monroe County in
Eastern Pennsylvania near the New Jersey border (Figure 1). The watershed drains into one
of the main tributaries of Delaware River and has very good water and biological quality, and
is designated as Special Protection Waters by the State and the Delaware River Basin
Commission (DRBC). Wild brown trout population, tourism and significant outdoor
recreation are major economic drivers of the area. Monroe County has the second fastest
growing population in the state of Pennsylvania. The county is threatened by high population
growth because of its attractive, pristine natural resources, and its proximity to the New York
City and Philadelphia metropolitan regions. Since 1980 the population of Monroe County has
nearly doubled and is projected to grow an additional 60% by 2020. Potential impacts,
among others, include a degradation and loss of the forested and agricultural lands and
deterioration of the local quality of life. Specifically, the concern is that the projected growth
and land use changes along with the accompanying increased groundwater withdrawals in the
watershed could well exceed sustainable levels, depleting groundwater and streamflows, and
resulting in the loss of the Creek's wild brown trout. Anticipated increase of imperviousness
due to urbanization in the watershed and projected increase in the demand for ground water
due to the anticipated population growth threaten the sustainability of current base flows.
Further, reduced infiltration and increased runoff rates have the potential to elevate peak
flows and increase flood hazards during large storm events.
The Delaware River Basin Commission in collaboration with the U.S. Environmental
Protection Agency (USEPA), Broadhead Watershed Association, Monroe County
Conservation District, Monroe County Planning Commission, U.S. Geological Survey
(USGS), and other stakeholders formed a consortium to study the potential future impacts of
the projected growth and land use changes on the sustainability of the natural resources in the
Pocono Creek watershed. The goals are to manage flows and growth in the Pocono Creek
watershed such that natural resources are sustainable. The goals will be partly achieved
through three integrated model studies by the USEPA, the U.S. Geological Survey (USGS),
and the PA Fish and Boat Commission (PA F&B). These model studies will evaluate the
effects of growth and land use change on ground water, streamflow, and the habitat of
Pocono Creek. This report documents the development of the Pocono Creek Watershed
model by the USEPA using the Soil and Water Assessment Tool (SWAT), its application to
assess the impact of projected urbanization on streamflow characteristics, and identification
of critical areas within the watershed having major contributions to changes in the
streamflow. The results of this model study will be linked to a USGS groundwater flow
model (MODFLOW) and the Pennsylvania Instream Flow Model (PIFM).
Distributed watershed models are utilized to better understand the role of hydrological
processes that govern surface and subsurface water movement. They provide tools for
environmental decision making and water resources planning and management. They are not
-------
only helpful in making model predictions of future flow conditions, but also in assessment of
hydrologic impacts of management measures scenarios, land cover and climate changes.
(961.614) B
B (961.615)
* Mount Pocono
climate station
(962,614)
(961,613
4 Kilometers (965,612)°
Streams
NEXRADcell centroicls
Elevation (meters)
| 1 183-219
219 - 254
254 - 290
290 - 326
326 - 362
362 - 397
| | 397 - 433
||433-469
| | 469 - 504
| | 504 - 540
^B 540 - 576
^B 576 -612
• 612 - 647
(966,614)
Q
*
Stroudsburg
climate station
Stroudsburg
USGS
stream gauge
Figure 1. Pocono Creek Watershed: Location, topography, USGS stream gauge, climate
stations and NEXRAD cell centroids (adapted from Kalin and Hantush, 2006a).
1.2 Objectives
The goal of the current model study and those by the USGS and the PA F&B, is to provide
the necessary science so that Monroe County and stakeholders can make informed decisions
that will assist them in developing and implementing sustainable water resources
management strategies without compromising the integrity of the watershed's habitat. To
-------
accomplish this overall goal, a watershed model is developed by the U.S. EPA to achieve the
following objectives:
1- Calibrate and validate a watershed model for Pocono Creek Watershed and examine
model performance with Next Generation Weather Radar (NEXRAD) data against that
based on surface raingauge precipitation data: The distributed hydrologic Soil Water
Assessment Tool (SWAT) will be calibrated and validated. Radar generated precipitation
estimates such as from NEXRAD products have found increasing usage in the hydrologic
community lately as an alternative source to gauge data. NEXRAD data can provide
information about the spatial distribution of precipitation patterns. If proven successful,
NEXRAD data can be a cost-effective alternative to the rather more costly and sparse
raingauge data.
2- Predict the impact of projected land use changes on annual average recharge distribution:
Spatial distribution of annual groundwater recharge rates will be computed pre- and after
urbanization build out for use by a USGS MODFLOW model to simulate the impact of
projected increase in groundwater withdrawals on base-flow reductions.
3- Predict the impact of projected land use changes on monthly median daily flows:
Monthly median daily flows will be computed pre and post urbanization build out for use
by the PIFM model which will be used to establish a relationship between flow
reductions and potential habitat degradation in the Pocono Creek watershed.
4- Evaluate model predictive uncertainty: It is a common practice to calibrate and validate
hydrologic and water quality models, but their forecasting abilities are rarely rigorously
evaluated. In this modeling effort, the extra step of evaluating model error propagation
will be conducted using time series analysis and Monte Carlo (MC) type simulations.
5- Compute the effect of urbanization on streamflow characteristics: It is anticipated that
projected population growth and urbanization in Pocono Creek will be accompanied by
an increase in impervious fraction of the watershed. The impact of land-use changes on
low, high, monthly average, and median flows will be particularly investigated, along
with 95% confidence band of the computed changes in flow characteristics.
6- Identify critical areas in the watershed: Once potential changes in flow characteristics are
predicted by the watershed model due to projected land-use changes, the inner
catchments contributing most of the changes will be identified using an index-based
methodology.
1.3 Report Organization
This report is organized as follows. In Section 2 the rationale for the selection of the SWAT
model to investigate the hydrology of Pocono Creek is established. A summary of model
representation of the most important components of the hydrologic cycle and related
parameters is presented. This summary is intended to define key model parameters and
provide insights into the process of model calibration. Section 3 describes the study area,
-------
model data and sources, and model calibration and validation using raingauge and NEXRAD
precipitation data. The viability of NEXRAD technology as an alternative to raingauge
networks as a source of spatio-temporal precipitation is established in the section. Error
propagation and model forecast capability is investigated in Section 4 using time series
analysis and Latin Hypercube Monte Carlo (MC) simulations. Section 5 explores potential
impacts of projected build out in the watershed on changes in key streamflow characteristics
using MC simulations. The impact of land use changes on (low, medium, and high) flow
frequency and duration is also investigated. Section 6 applies an index approach to identify
and rank subwatersheds or areas within subwatersheds that may contribute mostly to
predicted changes in streamflow characteristics. Each section ends with conclusions. The
final summary and conclusions of this model study appear in Section 7. References are
included in Section 8.
1.4 Summary
Pocono Creek watershed is threatened by high population growth and urbanization. Potential
impacts include degradation and loss of the forested and agricultural lands and deterioration
of the local quality of life. The condition of the wild brown trout habitat has been identified
as an indicator of the health of the watershed. Of specific concern is the potential impact of
projected population growth and land use change on the reduction of base flow and the
impact this may have on the degradation of wild brown trout habitats in Pocono Creek.
Projected increase in the imperviousness in the watershed, on the other hand, is expected to
increase flood frequency and reduce the recurrence interval of high flows.
Upon request by the DRBC and the EPA Region III, the National Risk Management
Research Laboratory of the USEPA has been tasked to investigate through a watershed
model study the impact of projected land use changes on potential alterations in the
hydrology of Pocono Creek. Constrained by the existing budget and available resources, it
was concluded that the SWAT model is a suitable and effective tool to conduct this modeling
study.
The objectives of the watershed model study are primarily three-fold. First, to calibrate
and validate a watershed model for Pocono Creek and compute the spatial groundwater
recharge distribution for use by a USGS groundwater flow model (MODFLOW) being
developed to assess the impact of land-use changes on base flow. Secondly, to compute
monthly median daily flows for use by the PIFM model that will be used to establish a
relationship between flow reductions and potential habitat degradation in Pocono Creek.
Thirdly, to investigate the impact of land-use changes on flow duration characteristics and
identify critical areas in the watershed. An important component of this study is investigating
NEXRAD as an alternative source of precipitation data, and the error propagation analysis
required to examine model forecast quality.
-------
2 Model Selection
2.1 Background
The distributed hydrologic Soil Water Assessment Tool (SWAT) (Neitsch et al., 2002a/b)
was chosen to fulfill the project objectives. The SWAT model was originally developed to
quantify the impact of land management practices in large, complex watersheds with varying
soils, land use, and management conditions over a long period of time, on the order of years.
It is developed and maintained by the US Department of Agriculture (USDA) scientists and
is freely available from . Although SWAT is mostly based
on empirical equations and simplified mass balance relationships, it is a widely-used model
and numerous applications can be found in the peer reviewed literature. For instance, as of
February 2006, the SWAT model web site cited 211 peer reviewed publications in the form
of journal papers or book chapters (). Borah (2002) reviewed eleven continuous-simulation and single-event
watershed scale models including SWAT. The study provides a better understanding of the
mathematical bases of the models. Kalin and Hantush (2003) reviewed key features and
capabilities of widely cited watershed scale hydrologic and water quality models, and
identified SWAT as one of the most suitable models for applications related to watershed
management. Robustness of SWAT for simulating watershed responses has also been
demonstrated in comparative studies by Saleh and Du (2002) and Van Liew et al. (2003).
Arnold and Fohrer (2005) provided a list of SWAT applications in the USA and worldwide.
The SWAT model development, operation, limitations, and assumptions were discussed
by Arnold et al. (1998). Srinivasan et al. (1998) reviewed the applications of the SWAT
model in streamflow prediction, sediment and nutrients transport, and effects of management
practices on water quality. Arnold and Allen (1996) evaluated the performance of different
hydrologic components of the SWAT model for three watersheds in Illinois (100-250 km2).
Comparing the model outputs to measured data, the calibrated model reasonably simulated
runoff, groundwater, and other components of hydrologic cycle for the study watersheds.
Most of the simulated average monthly outputs were within 5% of the historical data and
nearly all of them were within 25%. The coefficient of determination (R2) was used to assess
the correlation between the observed and simulated average monthly variables. Also, the
interaction among various components of hydrologic budgets was recognized to be realistic.
SWAT was utilized in a study by Arnold et al. (2000) to compare the performance of two
baseflow and groundwater recharge models. The first model was the water balance
components of the SWAT model. A combination of a digital hydrograph separation tool and
a modified hydrograph recession curve displacement technique composed the second model.
The results of the two models were in general agreement in the Upper Mississippi river basin.
A detailed procedure for calibration of SWAT was laid out by Santhi et al. (2001). Jha et al.
(2003) found curve number (CN) as the most sensitive parameter in streamflow prediction
with SWAT. Muleta and Nicklow (2005) applied SWAT coupled with automated calibration
to estimate daily flow and sediment yield in a 133 km2 Southern Illinois watershed. Eckhardt
and Arnold (2001) developed a version of SWAT having global optimization algorithm
(SWAT-G) to model daily flow in an 81 km2 watershed in Germany. Fohrer et al. (2002)
used SWAT-G in conjunction with two other GIS based agricultural economy and ecology
-------
models in a mountainous 60 km2 watershed in Germany to analyze the effect of land use
changes. Sophocleaous and Perkins (2000) integrated SWAT with MODFLOW and applied
the integrated modeling system to three different Kansas watersheds. Tripathi et al. (2004)
showed on a 92.5 km2 Indian watershed that SWAT can successfully simulate average annual
and monthly flow and sediment yield even if weather input is obtained through SWAT's
weather generator. One may also refer to Jayakrishnan et al. (2005) on SWAT applications to
water resources management. In conclusion, SWAT performance has been extensively
validated for streamflow, and sediment and nutrients yield predictions for different regions of
United States and outside.
2.2 Summary of Hydrologic Component of SWAT Model
SWAT is a distributed, deterministic process-based hydrologic model (Neitsch et al.,
2002a/b). The AVSWAT (Di Luzio et al., 2002) graphical user interface (GUI) which runs
under ArcView GIS is used to preprocess model data, run the SWAT model, and post
process model outputs. SWAT uses readily available inputs and has the capability of routing
runoff and chemicals through streams and reservoirs, adding flows and input measured data
from point sources, and is capable of simulating long periods for computing the effect of
management changes. Major components of the model include weather, surface runoff,
return flow, percolation, evapotranspiration (ET), transmission losses, pond & reservoir
storage, crop growth & irrigation, groundwater flow, reach routing, nutrient & pesticide
loading, and water transfer.
Input data needed to run the SWAT model includes soil, land use, weather, rainfall,
management conditions, stream network, and watershed configuration data. AVSWAT has
the capability of extracting most of these model parameters from readily available GIS maps
such as digital elevation models (DEM), land use maps, STATSGO soil maps, etc. Below is
a short summary of the SWAT model from the model manual and theoretical documentation
version 2000 (Neitsch et al., 2002a; 2002b).
SWAT partitions the watershed into subunits including subbasins, reach/main channel
segments, impoundments on main channel network, and point sources to set up a watershed.
Subbasins are divided into hydrologic response units (HRUs) which are portions of subbasins
with unique land use/management/soil attributes. AVSWAT enables extraction of input
parameters easily. It uses Digital Elevation Models (DEM) as input to extract the channel
network and delineate the watershed and subwatersheds, the resolution of which depends on
the user provided threshold area which is required to initiate a first order channel. The
threshold area can be chosen in such a way that the resultant channel network resembles the
one provided in topographic maps. The user needs to provide two threshold values to create
HRUs, one for land use and one for soil. Land uses that cover a percentage of the subbasin
area less than the threshold level are considered minor and thus eliminated. After the
elimination process, the areas of the remaining land uses are reapportioned so that 100% of
the land area in the subbasin is modeled. The soil threshold is applied next in a similar
fashion to eliminate minor soil types that occupy negligible portions of the HRUs.
In SWAT, the land phase of the hydrologic cycle is based on the water balance equation:
-------
SWt = SW0+^ (Rday,, - <2OTr/, - EaJ - Wseepj - <9W ) (2.1)
;=1
where SWt is the final soil water content (mm), SWo is the initial soil water content (mm), t is
the time (days), Rday,i is the amount of precipitation on day /' (mm), Qsurf,i is the amount of
surface runoff on day /' (mm), Ea,i is the amount of evapotranspiration on day / (mm), wseep,i is
the amount of percolation and bypass flow exiting the soil profile bottom on day /' (mm), and
Qgw,i is the amount of return flow on day / (mm).
Snowmelt is included with rainfall in the calculation of runoff and percolation; it is
controlled by the air and snow pack temperature, the melting rate, and the areal coverage of
snow. The mass balance for the snow pack is given by:
SNOt = SNOo + Rday - Esub - SNOmit (2.2)
where SNOt is the water content of the snow pack at the end of a day (mm), SNOo is the
initial water content of snow pack (mm), Rday is the amount of precipitation on a given day
(mm), Esub is the amount of sublimation on a given day (mm), and SNOmit is the amount of
snowmelt on a given day (mm). This equation assumes that the water released from
snowmelt is evenly distributed over the 24 hours of the day.
SWAT uses the SCS curve number method (USDA Soil Conservation Service, 1972) or
Green & Ampt infiltration method (Green and Ampt, 1911) to compute surface runoff
volume for each HRU. The former option is utilized in this study. The SCS runoff equation is
an empirical model that was developed to provide a consistent basis for estimating the
amounts of runoff under varying land use, soil types, and antecedent moisture conditions
(Rallison and Miller, 1981). The SCS curve number equation is:
(Rda-0.2S)2
day
where, S is the retention parameter (mm). In this equation the initial abstraction, which
includes surface storage, interception and infiltration prior to runoff, is approximated as 0.25".
The retention parameter is defined as:
CN
(2.4)
^ }
in which C7Vis the curve number for the day, which is a function of the soil's permeability,
land use and antecedent soil water conditions.
Evapotranspiration (ET) is the primary mechanism by which the water is removed from a
watershed. It includes all processes by which water at the earth's surface is converted to
water vapor: evaporation from the plant canopy, transpiration, sublimation and evaporation
from the soil. SWAT calculates actual ET from potential evapotranspiration (PET). The latter
is estimated by three methods in SWAT: the Penman-Monteith method, the Priestly-Taylor
-------
method, and the Hargreaves method. The Penman-Monteith method requires solar radiation,
air temperature, relative humidity and wind speed. The Priestly-Taylor method requires solar
radiation, air temperature and relative humidity. The Hargreaves method requires air
temperature only. Penman-Monteith is used in this study.
Once surface runoff is calculated, the amount of surface runoff released to the main
channel is computed from
—surlag
l-e 'conc
Qch,i ~ \Qsurf,i + Qstor,i-l)'
(2.5)
where Qch,t is the amount of surface runoff discharged to the main channel on day /' (mm),
Qsurf.t is the amount of surface runoff generated in the subbasin on day /' (mm), Qst0r,i-i is the
surface runoff stored or lagged from day /'-I (mm), surlag is the surface runoff lag
coefficient, and tconc is the time of concentration for the subbasin (hrs). The last expression
within the bracket on the right hand side of Equation (2.5) represents the fraction of total
available water allowed to enter the reach on a given day. Remaining water becomes
available water for the next day (Qstor.i).
The movement of water through the channel network of the watershed to the outlet is
routed in main channel and reservoirs. Flow is routed through the channel using a variable
storage coefficient method developed by Williams (1969) or the Muskingum routing method,
the latter of which is employed in this study. Transmission losses, which reduce runoff
volume as the flood wave travels downstream, are also accounted for by the model. SWAT is
a continuous time model, i.e., a long-term yield model. The model is not designed to
simulate detailed, single-event flood routing.
The water balance for the shallow aquifer is:
«&/,,/ = «&/,,/-! + Wrchrg,l ~ Qgw,i ~ Wrevap,, ~ H 'deep,, ~ Wpump,Sh,, (2-6)
in which aqsh,t is the amount of water stored in the shallow aquifer on day / (mm), aqsh,t-i is
the amount of water stored in the shallow aquifer on day /'-I (mm), wrchrg,i is the amount of
recharge entering the aquifer on day / (mm), Qgw>i is the groundwater discharge, or base flow,
into the main channel on day /' (mm), wrevap>i is the amount of water moving into the soil zone
in response to water deficiencies on day /' (mm), Wdeep,i is the amount of water percolating
from the shallow aquifer into the deep aquifer on day /' (mm), and wpump:Sh,i is the amount of
water removed from the shallow aquifer by pumping on day /' (mm). The recharge to the
aquifer on a given day is calculated:
Wrchrg,i = - e~sW WSeep,, + C~S" • Wrchrg,,-l (2-7)
where Sgw is the delay time or drainage time of the overlying geologic formations (days),
p,i is the total amount of water exiting the bottom of the soil profile on day /' (mm), and
i is the amount of recharge entering the aquifer on day /-I (mm).
-------
Base flow is computed by SWAT using this equation:
\+ \
(2.8)
where Qgw,i-i is the groundwater flow into the main channel on day /'-I (mm), agw is the base
flow recession constant, and At is the time step (1 day).
2.3 Summary
The rationale for the selection of SWAT was laid down in this section and its track record of
applications to watershed management was discussed. It should be noted that while it is not
the policy of the U.S. EPA to promote the use of a particular model, the unavailability of a
detailed survey of Pocono Creek channels' geometry and characteristics and lack of
information about bathymetric and hydraulic characteristics of impoundments, ponds, and
wetlands in the watershed precluded the use of more physically based, complex watershed
models that are currently available to the public. It will be evident throughout the analyses
hereafter that the selection of the SWAT model was indeed an appropriate decision.
SWAT is a distributed, process-based watershed model, but with significant number of
empirical relationships. It is one of the most suitable models for assessing the impact of
management practices and land disturbances on watershed responses, and has a solid track
record of applications.
-------
3 Watershed Model
3.1 Study Area
The 120 km2 Pocono Creek watershed is located between the latitudes 40°59'N-41°06'N and
longitudes 75°14'W-75°26'W in Monroe County, Eastern Pennsylvania near the New Jersey
state border (Figure 1), located within the Delaware River Basin. Pocono Creek's 26 km long
watershed valley drains from the Pocono Plateau in its headwaters eventually into the
Brodhead Creek, a tributary to the Delaware River. The model is constructed for the area
upstream from the USGS streamflow gauge station that is located about 6.4 km upstream
from the mouth near the city of Stroudsburg, PA (Figure 1) and drains an area about 98.87
km2.
Table 1 summarizes the data used in this study along with their sources and formats.
Precipitation, temperature, humidity and wind speed data were obtained from two National
Weather Service (NWS) climate stations; Mount Pocono to the North and Stroudsburg to the
East. As can be seen from Figure 1, both stations are outside the watershed boundary. To
study the potential effects of this on model performance, NEXRAD data were also utilized
during the study as an alternative precipitation data source. Specifically, the XMRG products
produced by The Middle Atlantic River Forecast Center (MARFC
) were used. XMRG precipitation files are generated in a
specific file format after analyses from both gauges and radar with some manual quality
control and are available at approximately 4 km cell resolutions. The small squares with dots
inside in Figure 1 represent the locations of the centroids of the NEXRAD precipitation cells.
The XMRG files for the MARFC region can be downloaded from
. Climate data from 1960 to
2004 indicate that, on average, annual precipitation is 1237 mm, and varies from a minimum
of 76 mm in February to a maximum of 125 mm in September. The temperature typically
varies from a minimum of-11 °C in January to a maximum of 26 °C in July.
The current land cover (Figure 2a) is dominantly forest (89%). Pasture constitutes about
3.5% and minor agricultural activities less than 0.2%. Residential, commercial and
transportation areas comprise about 5.8% of the watershed including the commercially
developed Route 611 corridor, Big Pocono State Park, Camelback Ski Area, the Nature
Conservancy's Tannersville Cranberry Bog, and state gamelands. Silt loam is the major soil
type in the watershed covering about 84.9%. Two other soil types in the watershed are sandy
loam (11.4%) and loam (3.7%). The elevation in the watershed changes from 183 m at the
outlet to 648 m near the Camel Back Ski Area. The topography in the watershed generally
has an average slope of 11% and ranges from 4% to 23%. A 30 m resolution digital elevation
model (DEM) is used in extraction of the stream network and delineation of the watershed
and sub watersheds. The watershed is divided into 29 subbasins. The STATSGO soil database
was used to acquire the soil-related model parameters (Figure 2b). Hydrologic response units
(FtRU) are generated by using 10% and 0% thresholds for land use and soil maps that
resulted in a total of 129 FIRUs.
10
-------
Table 1. Data Used in the Model Construct and Their Sources
Data
Source
Additional Info
Soil
Land Use
Elevation
Climate
Radar
Stream
flow
SWAT build in
DRBC
DRBC provided, also available at
http://www.pasda.psu.edu/
NOAA National Data Center
http://nndc.noaa.gov
NWS Hydrologic Data Systems
Group:
http://dipper.nws.noaa.gov/hdsb/data
/nexrad/nexrad. html
http://waterdata.usgs.gov/pa/nwis/uv
?dd_cd=02&format=pre&period=16
&site no=01441495
USDA STATSGO soil database
National Land Cover Data (NLCD)
30m resolution DEM
Stations: Mount Pocono, 4F08'N /
75°23'W; Stroudsburg, 4F01'N/
75°11'W; hourly (daily for Stroudsburg)
precipitation, temperature, wind speed,
relative humidity
NEXRAD, Multisensor Precipitation
Estimator (MPE) Data in XMRG format.
Hourly precipitation at approximately 4
km resolution.
USGS 01441495 Pocono Creek ab
Wigwam Run near Stroudsburg, PA, Lat
40°59'27", long 75°15'20", data collected
every 15 min.
3.2 Base flow Separation
Stream flow is usually partitioned into two parts: the fast and the slow response components,
the latter of which is due to the base flow contribution. Any other contribution to stream flow
by various mechanisms can be deemed as the fast response component. SWAT computes the
base flow and surface runoff components of the stream flow separately. Although some
parameters play an interactive role such as the CN, some parameters only affect one
individual component of the flow. For instance, Manning's roughness only affects surface
runoff, whereas base flow recession constant only affects base flow. Hence, to better
calibrate model parameters it is necessary to partition the observed stream flow hydrograph
into base flow and surface runoff components, both of which then become estimated
quantities rather than observed.
We estimated base flow using the method described in Arnold et al. (1995) and Arnold
and Allen (1999). This is a recursive digital filter technique which was originally used in
signal processing and filtering. The filter can be passed over the stream flow data several
times with each pass resulting in less base flow. Various other methods of base flow
separation are also available and each can give significantly different base-flow estimates
which clearly affect not only model calibrated parameters but also model performance.
11
-------
(a)
Open Water
Perennial Ice/Snow
Low Intensity Residential
High Intensity Residential
Commercial/Industrial
Bare Rock/Sand/Clay
Quarries/Strip Mines
| | Transitional
Deciduous Forest
Evergreen Forest
Mixed Forest
Shrubland
| | Orchards/Vineyards
Grasslands/Herbaceous
Pasture/Hay
Row Crops
Small Grains
Fallow
Urban/Recreational Grasses
Woody Wetlands
Emergent Herbaceous Wetlands
(b)
Statsgo soils
PA024
PA025
PA026
PA027
PA030
PA031
PA037
Figure 2. Pocono Creek watershed: a) Land use map; b) Soil map (STATSGO).
3.3 Next Generation Weather Radar (NEXRAD)
Both the Mount Pocono and Stroudsburg climate stations are located outside the watershed
boundary (Figure 1) and provide point measurements. This raises the question whether the
data from these two gauge stations are representative of spatial rainfall pattern over the
Pocono Creek watershed which is partly covered with mountains. It is well known that
precipitation may change significantly in mountainous areas due to orographic effects.
Precipitation is usually higher in upper elevations. To be able to realistically answer this
question several raingauges inside the watershed boundary are needed so that spatial
distribution of the precipitation pattern can be attained. The Next Generation Weather Radar
(NEXRAD) precipitation data offer an opportunity to overlay spatially distributed
precipitation over the entire watershed. NEXRAD also provides a finer temporal resolution
(hourly) compared to daily values from the climate stations that are available. In any case,
hourly precipitation is not required for this study and daily precipitation suffices for the
current model application. Because the ultimate goal is modeling base flow and stream flow
in the Pocono Creek watershed, gauge driven model performance is compared to NEXRAD
driven ones. This sheds light on whether or not consideration of the distributed nature of
precipitation improves the model performance, and whether NEXRAD can be relied on as an
alternative to the surface raingauge measurements. For further details on processing the
NEXRAD data for use in a watershed model, the reader may refer to Kalin and Hantush
(2006a).
3.4 Model Calibration
Models are only simplified representations of natural processes. Even fully physically-based
models cannot avoid simplifications. In addition, parametric uncertainty and measurement
errors make model calibration inevitable in most modeling exercises. A split data set
approach is implemented in this study to calibrate and validate the SWAT model. The period
12
-------
from 7/1/02 to 5/31/04 of the daily flow data, aggregated into monthly flows, is used for
calibration and the remaining data from 6/1/04 to 4/30/05 is used for validation. Three to five
years of data are typically required in calibration, although along with a good set of data and
proper objective function, a single year of data has been shown to be adequate (Sorooshian et
al., 1983). The statistical measures of Mass Balance Error (MBE), coefficient of
determination (R2) and Nash-Sutcliffe (1970) efficiency (ENS) are used as indicators of model
performance that are defined as
N N
YO -Vo, . — -
/_^z^sim,i /_^z^oos,i f^. — r)
i=\ ;=1 _ xi sim z--o
JL n
Vn ^obs
(3.1)
R2=^
l2
(3.2)
Z^ (OSim,i
i=\
(o o\
3.3)
•th
where Qsim,i and Q0bs,i are simulated and observed or estimated flows at i observation,
respectively, N is the number of observations. Similarly, Osim and Oobs are the average
simulated and observed flows over the simulation period. The coefficient of determination
describes the proportion of the total variances in the observed data that can be explained by
the model and ranges from 0 to 1, whereas Nash-Sutcliffe efficiency is a measure of how
well the plot of observed versus predicted values fit the 1:1 line, and can vary from -GO to 1.
A negative ENS indicates that model predictions are not better than the average of observed
data.
The calibration process is adapted from Santhi et al. (2001) and is summarized in Figure
3. The threshold values set for MBE, R2 and ENS in the figure are at monthly time scales.
Curve numbers (C7V) of each soil/land use combination are calibrated first to meet the criteria
set for surface runoff (SR). At the next stage the GW REVAP coefficient which is a limiting
factor for the maximum amount of water that can be removed from the aquifer to the
overlying unsaturated zone due to moisture deficit, threshold depth of water in shallow
aquifer required for base flow to occur (GWQMN), delay time for aquifer recharge
(GW_DELAY), and soil evaporation compensation factor (ESCO) are calibrated to meet the
stream flow (SF) and base flow (BF) criteria. Once the model is calibrated using measured
stream flows for the calibration period, the model is validated using the data for the above
stated validation period and MBE, R2 and ENS are checked against the threshold values shown
in Figure 3. The threshold values set for MBE, R2 and ENS statistical measures must be
13
-------
satisfied for a successfully calibrated model. Before calibrating the model parameters, an
automated sensitivity analysis was performed with a version of SWAT called AVSWATX
that is based on Latin Hypercube (LH) and One factor At a Time (OAT) sampling. The
sensitivity analysis revealed that CN is the most sensitive model parameter (Kalin and
Hantush, 2006a).
Separate SR and BF for
measured daily flow
If avrg. of sim. SR is within ± 15%
of avrg. measured SR and
R2>0.6, EN_S>0.5
YES
If avrg. of sim. SF is within ± 15%
of avrg. measured SF and
R2>0.6, EN_S>0.5
NO
Adjust GW_REVAP,
GWQMN, ESCO,
GW DELAY
YES
If avrg. of sim. BF is within
± 15% of average measured BF
NO
Adjust GW_REVAP,
GWQMN, ESCO,
GW DELAY
YES
Calibration complete
•SR: Surface Runoff
•SF: Stream Flow
•BF: Base Flow
•R2: Coefficient of determination
[0,1]
•EN.S: Coefficient of efficiency or
Nash-Sutcliffe Efficiency [-°°,1]
•CN: Curve Number
•GW_REVAP: Groundwater
revap* coefficient
•GWQMN: Threshold depth of
water in shallow aquifer required
for return flow to occur
•ESCO: Soil Evaporation
compensation factor
•GW_DELAY: Delay time for
aquifer recharge (days)
* movement of water from
unconfined aquifer to overlying
unsaturated layers
Figure 3. Calibration procedure in the SWAT model (adapted from Santhi et al., 2001).
Figure 4 compares the daily and monthly areally averaged NEXRAD estimated
precipitation to measurements at the gauges. All the values in the figure are in mm. Mount
Pocono station falls in between two NEXRAD cells (Figure 1): (961,614) and (961,615)
where the numbers in parentheses represent the X and Y coordinates of the centroids of
NEXRAD cells in the FtRAP coordinate system. Therefore, we compared Mount Pocono
point precipitation measurements to NEXRAD estimates at both cells as well as to their
arithmetic averages. The Stroudsburg station falls inside the grid (966,614) as seen in Figure
1. Estimated hourly precipitations were
14
-------
150
CC
I
400
5 300 -
SJ-
= 200 -
0)
100 -
Monthly-Scale
(1/1/02-4/30/05)
y = 1.03x + 8.96
R2 = 0.87
100 200 300
Mt. Pocono Gauge
150
400
150 -,
Daily Scale
120 -\ (1/1/02-4/30/05)
CD
82-
90 -
60
30
y=1.00x+0.25
R2 = 0.91
150
30 60 90 120 150
Mt. Pocono Gauge
30 60 90 120
Mt. Pocono Gauge
400
in
5 300
5
S-
1 20°
o
CL
2 100
I
Monthly-Scale
(1/1/02-4/30/05)
y=1.01x + 5.21
R2 = 0.94
100 200 300
Mt. Pocono Gauge
350
400
150
30 60 90 120
Mt. Pocono Gauge
150
Q_ '
<
or .
400
300 -
in"
^- 200 -
5
2L
100 -
0
Monthly-Scale
(1/1/02-4/30/05)
y=1.02x+7.09
R2 = 0.92
0 100 200 300 400
Mt. Pocono Gauge
^300-
s- 25° -
CD
a 200 -
8 150 -
< 100-
1 50-
0 -
Monthly-Scale
(1/1/02-2/28/05)
y=0.95x + 1.87
R2 = 0.90
30 60 90 120
Stroudsburg Gauge
150
50 100 150 200 250 300 350
Stroudsburg Gauge
Figure 4. Comparison of NEXRAD estimated precipitation with measurements from local
gauge stations. Values shown are in mm.
aggregated in these cells to obtain the daily and monthly precipitation estimates. We
computed daily values from 7:00 AM to 7:00 AM at the cell (966,614) to be consistent with
the reported daily precipitation data at the Stroudsburg gauge station. Figure 4 reveals that
(961,615) represents Mount Pocono relatively better than both (961,614) and the average of
the two, as it has a higher R2, the slope of its regression equation is closer to 1 and has a
smaller intercept. Further, (961,615) overestimates precipitation by only 5.7% compared to
10.8% overestimation of (961,614) from 1/1/2002 to 4/30/2005. Over the time period from
1/1/2002 to 2/28/2005, NEXRAD underestimates precipitation in the Stroudsburg station by
10.0%. In this watershed, comparisons with gauge measurements indicate that NEXRAD
technology provides an alternative source of precipitation data. Note that NEXRAD
estimates are areal average precipitations of an approximately 4x4 km2 grid in contrast to
point measurements (-100 cm2) of gauge stations. Differences between the two are partly
attributed to variations of point measurements from areal average estimates of precipitations,
and partly due to measurement errors associated with the instrument itself.
15
-------
SWAT is calibrated both with rain gauge and NEXRAD as the precipitation data source.
The simulations are performed starting from 7/1/1970, i.e., with a 32 years of warm up
period to minimize the effect of initial conditions. For more details on calibration using
raingauge data, the readers are referred to Kalin and Hantush (2006a). The surface runoff lag
coefficient (surlag) had to be adjusted to improve the daily simulation performances. The CN
values calibrated based on raingauge data are reduced in all HRUs by 1.5 (for this specific
site) when NEXRAD was the data source. The volume of precipitation over the study
watershed estimated with NEXRAD is greater than that estimated from raingauges, and thus
this adjustment in CN had to be made. Experimentation with other parameters revealed no
further improvement in the model performance. Figure 5 compares SWAT simulations to
observed SF and estimated SR at the monthly time scale and SF at the daily time scale only
for NEXRAD. Rain gauge simulations not shown on the figure as hydrographs are very
similar to NEXRAD generated hydrographs.
in
I 6
ce
0)
o
>2 o
i_ z
w
• Observed
- SWAT - nexrad
CD r- 00 O5 O
12
_(0 g
?
w
- Observed
-SWAT- nexrad
16
-------
40
-Observed
-SWAT-Nexrad
55
Figure 5. Measured stream flow (SF) and estimated surface runoff (SR) plotted against
SWAT simulated counterparts for the calibration (7/1/2002-5/31/2004) period.
NEXRAD estimates are used as precipitation data source. Top two panels show
monthly results, and the bottom panel shows daily results.
As Figure 5 shows, overall the model performs well at both time scales. These results are
comparable to gauge driven simulations (Kalin and Hantush, 2006a). The NEXRAD driven
model performance statistics MBE, R2 and ENS are summarized in Table 2 and are well within
the calibration threshold values shown in Figure 3. From the results obtained, although with a
relatively smaller calibrated CN, it is reasonable to conclude that NEXRAD estimated
precipitation data is a good alternative to gauge measured precipitation data. The benefit of
using distributed rainfall data was negligible at both time scales for the particular application.
On the other hand, noticeable improvement is expected for stream flow estimates at the
interior subwatershed outlets with the use of distributed rainfall data over the available
raingauge data. In the next section, an 11 month period split sample data set is used to
validate the calibrated model.
Table 2. Mass Balance Error (MBE), Coefficient of Determination, R2, and Nash-Sutcliffe
Efficiency, ENS for Stream Flow (SF), Base flow (BF) and Surface Runoff (SR) During the
Calibration Period with NEXRAD as Precipitation Data Source. Values in Parenthesis
Indicate Values when Rain Gauges are Used as Precipitation Data Source
Stream
MBE
Monthly -3.8
(-3.8)
Daily
Flow (SF)
R ENS
0.85 0.84
(0.85) (0.83)
0.74 0.73
(0.74) (0.74)
Base Flow (BF)
MBE
% R ENS
-4.5 0.31 0.05
(-4.0) (0.30) (0.08)
Surface Runoff (SR)
MBE
% R ENS
-3.2 0.79 0.79
(-3.6) (0.77) (0.77)
17
-------
3.5 Model Validation
The period from 6/1/2004 to 4/30/2005 is chosen as the validation period. Simulations are
performed again starting from 7/1/1970 to curtail the effect of initial conditions. Figure 6
compares observed values to model simulations conducted with NEXRAD data. In the figure
the top two panels compare monthly SF and SR and the bottom figure is for daily SF results.
From visual inspection, one can conclude that model simulations match well with observed
SF and estimated SR.
12
b 6
ce
0)
o
- Observed
-SWAT- nexrad
15
JO
CO
_§
1
LL
E
ro
0)
12 -
9-
6-
3
• Observed
-SWAT-nexrad
SSSSSSS8888
Figure 6. Measured streamflow (SF) and estimated surface runoff (SR) plotted against
SWAT simulated counterparts for the validation period (6/1/2004-4/30/2005).
NEXRAD estimates are used as precipitation data source. Top and bottom panels
denote monthly and daily results, respectively.
Table 3 summarizes the model efficiencies at the monthly and daily time scales, respectively,
validated with NEXRAD and raingauge data. In the table "g" denotes gauge input and "n"
indicates NEXRAD input precipitation. In all cases, the model calibration criteria are
satisfied. As can be seen, MBE in SR decreased in absolute value from 15.1% to 7.9% with
the NEXRAD input simulations. In SF and BF, MBE dropped in absolute value to 5.6% from
13.5% and to 1.2% from 10.6%, respectively, all underestimations. Overall NEXRAD input
18
-------
monthly simulations outperformed raingauge based simulations. Daily statistics showed
mixed results.
Table 3. Mass Balance Error (MBE), Coefficient of Determination, R , and Nash-Sutcliffe
Efficiency, ENS for Stream Flow (SF), Base flow (BF) and Surface Runoff (SR)
During the Validation Period, g and n Represent Gauge and NEXRAD Input
Precipitations, Respectively
Precip.
source
3 §
o
S n
g
'ci
C) n
Stream Flow
MBE
/O/ \ r>2
\/0) /v
-13.5 0.81
-5.6 0.89
0.70
0.66
(SF)
ENS
0.66
0.75
0.64
0.62
Base Flow (BF)
MBE
% R ENS
-10.6 0.13 -0.26
-1.2 0.06 -0.40
Surface Runoff (SR)
MBE
O/ 1)2 77
/o K &NS
-15.1 0.83 0.73
-7.9 0.84 0.77
It is interesting to note that while at the monthly time scale NEXRAD had better
performance statistics, at the daily time scale gauge driven simulations had slightly better
statistics. Several reasons may have contributed to this. First, calibration efforts were more
focused to monthly time scale as part of the project goal. Secondly, NEXRAD calibration
was built on the calibrated parameters with gauge driven data, therefore there is a small bias.
An additional set of 5 months streamflow data in the period 5/1/2005 to 9/30/2005, which
was not available during model validation, was downloaded from the source and was utilized
to further evaluate the model performance. The model was run with NEXRAD data and
corresponding calibrated parameters during this post-validation period, and the results are
plotted in Figure 7. Interestingly, expanding the validation period 6/1/2004 - 4/30/2005 to
9/30/2005 results in slightly improved results; the MBE in absolute value decreased from
2.9% to 0.8%, daily R1
0.64.
increased to 0.68 from 0.66, and daily ENS increased from 0.62 to
19
-------
Figure 7. Measured daily streamflows plotted against SWAT simulated counterparts for the
post-validation period (5/1/2005-9/30/2005). NEXRAD estimates are used as
precipitation data source.
3.6 Conclusions
In this section the hydrology in the Pocono Creek watershed was modeled using a distributed
parameter watershed model, and in particular the potential for using NEXRAD as an
alternative source of precipitation data to raingauge stations was explored. The SWAT model
was calibrated and validated for the Pocono Creek watershed.
The model was first calibrated using precipitation data from two gauge stations located
outside the watershed boundary. Model performance was evaluated with computed model
efficiency statistics, i.e., mass balance error (MBE), coefficient of determination (R2), and
Nash-Sutcliffe Efficiency (ENS). Simulation results were promising at the monthly and daily
time scales.
As an alternative data source, the use of radar generated precipitation data (NEXRAD)
was appraised. NEXRAD data obtained from the Middle Atlantic River Forecast Center
(MARFC) provided hourly precipitation estimates over approximately 4x4 km2 grids.
Measured precipitation values at the two gauge stations were in close agreement with the
NEXRAD estimated values.
The SWAT model was fed with the spatially distributed NEXRAD precipitation data and
recalibrated by reducing the average curve numbers (C7V) by a value of 1.5 that were
obtained from model calibration when gauge precipitation data was used as input. The
streamflow and surface runoff hydrographs were in close agreement with measured
hydrographs. Efficiency statistics MBE, R2 and ENS were close to ones obtained with gauge
driven simulations. The model was then run for an additional 11 month period for validation
purposes both with gauge and NEXRAD data. NEXRAD generated smaller MBE. The values
of R2 and ENs were higher at the monthly time scale with NEXRAD, whereas at the daily
time scale gauge driven simulations resulted in improved measures of fit. An additional 5
months set of streamflow data was used to further evaluate the performance of the already
validated model.
20
-------
This section showed that spatially distributed precipitation data obtained through radar
reflectivity measurements provide a viable alternative to raingauge measurements. However,
estimation of precipitation with the help of radar still needs improvements. At present, data
from raingauges will continue to be relied upon to correct NEXRAD hourly digital
precipitation for mean field bias. Future refinement of this technology may provide a cost-
effective alternative source of precipitation data, and may reduce the need for the costlier
raingauges for large scale watershed applications.
21
-------
4 Model Prediction Uncertainty
4.1 Introduction
It is a common practice that the sequence of sensitivity analysis —» calibration —» validation
is followed during applications of distributed hydrologic/water quality models, yet rigorous
attempts are rarely made to assess model predictive uncertainty. Once the model is calibrated
with one set of data and validated with another set, it is used for predicting the impact of
management practices, land-use changes, and/or long-term climate changes, however, with
much less regard to uncertainty bands of predictions. The process of examining model
forecast reliability (post-validation) is an important consideration in development of
watershed management plans. Modeling uncertainty should be rigorously addressed in
development and application of models, especially when the model outcome might have
implications on policy, watershed planning and management, and when stakeholders are
affected by the decisions contingent upon model-supported analyses (NRC, 2001).
Implications of model uncertainty should be factored in the decision making process.
The analysis of uncertainty associated with utility of simulation models appears mostly in
the scientific, research literature (e.g., Spear and Hornberger, 1980; Beven and Binely, 1992;
Spear et al., 1994; van der Perk and Bierkens, 1997; Saltelli et al., 2000; Hossain et al., 2004;
Carpenter and Georgakakos, 2004; and Pebesma et al., 2005). While various approaches exist
for estimating distributed watershed/water quality model prediction uncertainty, careful
examination of the scientific literature reveals three dominant approaches. Amongst the three
approaches, the Bayesian approach is perhaps the most poular. Examples include Bayesian
Monte Carlo (e.g., Dilks et al., 1992; and van der Perk and Bierkens, 1997) and Generalized
Likelihood Uncertainty Estimation (GLUE) (Beven and Binely, 1992; Freer et al., 1996;
Beven and Freer, 2001; Hossain et al., 2004; and Hossain and Anagnostou, 2005)
methodologies. A second approach relies directly on Monte Carlo (MC) method to obtain an
ensemble of model outputs by independently sampling model parameters from prior
distributions based on ranges (i.e., with minima and maxima) derived from literature, or
based on new information gained from experience and model calibration (e.g., Binley et al.,
1991; Carpenter and Georgakakos, 2004; and Hantush and Kalin, 2005). A third approach
can be discerned in which the model noise or residual error is accounted for explicitly (e.g.,
Sorooshian and Dracup, 1980; and Pebesma et al., 2005). In this approach, model predictions
and observations are given as time series data, and attempts are made to fit a deterministic or
statistical relationship to the residual time series. Among these three approaches, the third
one has been given the least attention. Although GLUE methodology has gained a wider
acceptance, the selection of a suitable relationship for computing relative likelihoods
associated with the ensemble simulations is a matter of choice, and therefore subjective in
that sense. Both the GLUE methodology and the third approach implicitly account for the
four major sources of uncertainty: model structural uncertainty, parametric uncertainty,
measurement uncertainty, and rainfall variability.
In this section, we fit a time series model to the residuals (observed - predicted) of daily
SWAT model output. The focus on model performance at the daily time scale is because
22
-------
predictions of monthly median daily flows are required inputs to the PIFM model for wild
brown trout habitats. The time series model is combined with Monte Carlo-type simulation
(Latin Hypercube) to estimate uncertainty bounds for model predictions. Using a new set of
streamflow data (i.e., different from calibration and validation data), SWAT model forecast
performance is evaluated by comparison of the 90% confidence interval or (uncertainty band)
with observed values. One of the advantages of the time series approach, which is described
in the following section, is that the commonly adopted Gaussian and statistical independence
assumptions for model errors are relaxed by brute force application of time series analysis
and using a nonparametric probabilistic approach. Further, MC-type simulation is conducted
on rather a simple time series model for the residual errors as opposed to the computationally
demanding watershed model.
4.2 Structural and Model Parameters Uncertainty
4.2.1 Prediction Error
Figure 8(a) shows the time series of model prediction errors (observed-simulated) based on
simulations during the calibration and validation period (7/1/2002-4/30/2005). The error is
defined by this equation
st=0t~Pt (
where &t is model output noise or prediction error; ot is measured streamflow; andpt is model
computed streamflow; and t is time index in days. The st accounts for model structural
uncertainty, parametric uncertainty, measurement uncertainty, and errors in the rainfall input.
Model structural uncertainty is usually associated with imperfect knowledge, and generally
referred to as epistemic uncertainty. A close look at the daily errors time series (Figure 8a)
reveals systematic errors shown as large positive or negative spikes that are not random in
nature. In many instances a large (-) error is immediately followed by a big (+) error. One
reason might be timing errors in recording either the precipitation or the streamflow.
Another, probably the more reasonable, explanation is the inability of the SWAT model to
use sub daily rainfall data. When a big rainfall event happens at the very end of a day, the
actual watershed response will be about 1 day delayed compared to what SWAT actually
simulates with daily rainfall data. This will clearly result in overprediction (negative error)
and underprediction (positive error) on two successive days, respectively. Hence, to
minimize this type of error to some extent, we decided to use three days (3-day) moving
average of the errors in constructing our model. While one may argue as to why daily
simulations are of interest, especially when SWAT is best suited for long-term simulations
(monthly and annual time steps), there are two compelling reasons for this consideration.
First, monthly median daily flows are a required input to the PIFM habitat model. Secondly,
there is not sufficient monthly streamflow data to conduct a meaningful time series analysis
to model errors. Therefore, in addition to smoothing the effect of the large errors during
significant events, a 3-day averaging of streamflows is intended to provide reasonable
estimates of daily flows for the wild brown trout habitat model. Further, sufficient data of 3-
day moving averages will be available to construct a time series model of residual errors et.
23
-------
Henceforth, st refers to model residual errors based on 3-day averages of ot andpt, and thus
deemed (i.e., st) a surrogate to (rather than exact) daily model prediction errors.
(a)
i/r 20
CO
I 10
0 -
co
T3
CD
"o
T3
£
CD
CD
O -20
CM
O
CM
O
CO
O
CO
o
CO
o
CD
CD
CD
5)
CM
CM
CM
O
in
o
CD
CM
in
CD
CD
CO
0
QJ O
o *=
II
9- CD
T3 >
CD CO
e
CD
(b)
20
10 -
0 -
O)
£ -10 -
-20
CM
O
CM
O
CO
CO
o
CO
o
CO
o
CNi
CO
o
O
CD
O
55
o
CNi
CM
i^
CM
CD
in
o
3
CM
in
o
o
CO
Figure 8. Time series of model prediction errors (Observed - SWAT simulated) in the
calibration and validation time period (7/1/2002-4/30/2005): (a) daily simulations,
(b) 3-day moving average.
Figure 8(b) shows the difference between 3-day average of observed minus predicted
times series of stream flows during the period (7/1/2002-4/30/2005). The errors are not
totally eliminated but their magnitudes are significantly reduced. The Minitab® statistical
computer package is used to identify a time series model and construct reasonable forecasts
for the future of st. In general, a multiplicative seasonal autoregressive integrated moving
average model denoted by AREVIA (p,d,q) x (PJ),Q)s is fitted to the series of interest, which
is the model prediction error st in our case. The parameters/*, d, and q are, respectively, the
24
-------
orders of the autoregressive, difference, and moving average operators; P, D, and Q are,
respectively, the orders of the seasonal autoregressive, seasonal moving average, and
seasonal difference operators; and S is the seasonal period. Autocorrelation function (ACF)
and partial autocorrelation function (PACF) are utilized and the procedure outlined by
Shumway (1988), which derive its basis from Box and Jenkins (1970), is followed to identify
the best model. The objective of the model identification process is to produce identically
distributed, independent residuals wt, with a minimum variance aw2 arising from fitting some
ARIMA (p,d,q) x (PJ),Q)s to the et time series.
The upper two panels in Figure 9 (i.e., Figure 9(a)) show the ACF and PACF of model
errors, et. It is obvious that the series display nonstationarity accentuated by the slowly
decaying ACF as a function of lag and the large positive value of PACF at lag 1. The ACF of
the first difference in Figure 9(b) contains a fairly strong negative peak (-0.44) at lag-3 and
zero thereafter, indicating that a seasonal moving average with S = 3 (days) and Q =1 might
be appropriate. The decreasing peaks of the PACF at multiples of 3 are due to the seasonal
moving average component. Therefore, ARIMA (0,1,0) x (0,0,1)3 appears to be, at the
moment, a suitable selection among, perhaps, other competing models. Table 4 lists several
ARIMA models with computed variance of residuals, cr^,, and various goodness-of-fit
measures: final prediction error (FPF), Akaike's information criterion (AIC), and Bayesian
information criterion (BIC). It is clear that the most appropriate model choice is ARIMA
(1,1,1) x (0,0,1)3, which has the smallest values of cP, FPF, AIC, and BIC. The model was
fitted to the error time series depicted in Figure 8(b) (recall the 3 -day averages of the actual
model errors), and the equation for determining the residuals wt is
-B)Vst = l-&B3(l-0 B)wt (4.2)
where Vst = st - st-i is the first difference operator; Bst = st-i is the backward shift operator;
& 9, and 0 are the fitted parameters. Equation (4.2) can be expanded to yield the following
representation for st:
et =l + ^-i -st-2 +wt -0w,_, -0w,_3 +00w,_4 (4-3)
25
-------
(a)
Autocorrelation Function for (Qbs-Sirn)
(with 5% significance limits for the autocorrelations)
Partial Autocorrelation Function for (Obs-Sim)
(with 5% significance limits for the partial autocorrelations)
1.0-
0.8-
0.6-
g 0.4-
o 0.2-
0
8 -n?-
3
< -0.4-
-0.6-
-0.8-
-1.0-
jjjhl1r
-- ^'JJlUJi1.'
1,0-
0,8-
c 0.6-
1 0.4-
t 0.2-
o
< -0.2-
~m
33 -0,4-
a -0.6-
-0.8-
-1.0-
|__i._4H__r_,.______l__r_. ___________________ --t--T-v
15 10 15 20 25 30 35 40 45 50 55 60 65 70 75 15 10 15 20 25 30 35 40 45 50 55 60 65 70 75
Lag Lag
(b)
Autocorrelation Function for 1st difference (Obs-Sim) Partial Autocorrelation Function for 1st difference (Obs-Sim)
(with 5% significance limits for the autocorrelations) (with 5% significance limits for the partial autocorrelations)
1.0-
0.8-
0.6-
§ 0.4-
0 0,2-
b 0,0-
a -0.2 -
< -0.4-
-0,6-
-0.8-
-1.0-
1,0
0,8
c 0.6
'iii 0.4
0
t 0.2
0 0.0
< -0.2
•S -0,4
S -0,6
-0.8
-1.0
1 '
15 10 15 20 25 30 35 40 45 50 55 60 65 70 75 15 10 15 20 25 30 35 40 45 50 55 60 65 70 75
Lag Lag
(c)
ACF of Residuals for (Obs-Sim) - (l,l,l)x(0,0,l)3 PACF of Residuals for (Obs-Sim) - (l,l,l)x(QAl)3
(with 5% significance limits for the autocorrelations) (with 5% significance limits for the partial autocorrelations)
0,2
.9
1 0-0
I
-0,2
r; if. ,. . 1 ... ,i , ........ r.i.
• ' |l | ' M "If "i| [' >'!'!' Ii"
r
i 0.2-
•J3
B
E
<
:"."il".""":.h,".":,"":"t""""T""".n":;.";""t"rr".
• 1 L • M •• l|. .-|| M II. II"
2
£ -0.2-
16 12 18 24 30 36 42 48 54 60 66 72
Lag
16 12 18 24 30 36 42 48 54 60 66 72
Lag
Figure 9. Autocorrelation (ACF) and partial autocorrelation (PACF) functions of: (a)
prediction error st, (b) first difference of residual error Vst = st - st.i, (c) residuals
(wt) for ARIMA( 1,1,1) x (0,0,1)3. The lag is in units of days.
26
-------
where
-------
to the past by preserving and reproducing the statistical characteristics of observed model
errors.
4.2.2 Nonparametric Probability Distribution
Nonparametric methods are applied when parametric probability distributions fail to describe
the stochastic sequence under consideration. Lall (1995) reviewed applications of
nonparametric function estimation in hydrology. A probability density function of arbitrary
shape can be locally approximated by a nonparametric model. This is particularly important
when commonly used probability distributions (e.g., normal, log-normal, exponential, etc.)
poorly fit the frequency of the observed stochastic series. However, despite their successful
applications to hydrology, nonparametric approaches have the limitation that they do not
extrapolate beyond the range of the record, since the sequence of future wtis hypothesized to
have a similar nonparametric functional form of the fitted probability density function (PDF).
In parametric methods, the density function is estimated by assuming that data are drawn
from a known parametric family of distributions. The methods of moments, maximum
likelihood estimation, or any other methods are commonly used to estimate the parameters of
the chosen PDF. In nonparametric approaches, a kernel function is often used to generalize
the density function estimation. Given a set of n observations w\, w2,..., wn, a mathematical
expression of a univariate kernel probability density estimator is ((e.g., Kim and Valdes,
2005)
(44)
where x is a random variable which stands for vc; K is a kernel function; n is number of
observations; and h is a bandwidth that controls the variance of the kernel function. Kim and
Valdes (2005) provide a list of kernel functions typically used in hydrology, the most widely
used one being the Gaussian:
(4.5)
An optimal estimate for the bandwidth for a Gaussian kernel is provided by Kim and Valdes
(2005) citing Silverman (1986),
(4.6)
Where cris the standard deviation of the observed record, which in this case is equal to the
standard deviation of observed residuals,
-------
Note that the residual term wt in the error stochastic model (4.3) is not directly observed in
the data sample sj, 82, ..., sn. This value can be approximated by assuming that w2 = w\ = w0
= w.i = 0 and then computing from (4.3)
(4.7)
for k = 3, 4, ..., n. Recall that the sequence si, S2, ..., sn, which is used to construct the time
series model (4.3), is computed from Eq. (4.1) using SWAT model simulations and
corresponding streamflow measurements during the calibration and validation period.
Although the residuals wt have zero-mean and are independent, they do not follow any of
the known parametric distributions. Figure 10 shows a relatively poor fit to the frequency
histogram by Gaussian PDF assuming wt~ N(0, cP). The nonparametric fit (4.4) provides a
remarkably improved local fit to the observed histogram.
1.5
JL
f(x)
N(x)
0.5
-3
-2
-1
0
birij,x
1
Figure 10. Frequency histogram, Gaussian PDF, and the nonparametric model fit to the
computed series of residuals Wk. Both the histogram of residuals and normal
distribution has zero mean and standard deviationov,.
In the following section we now combine nonparametric random generation and Latin
Hypercube Sampling (LHS) with the stochastic model (4.3) to synthesize the ensemble of
model error time series wt.
4.2.3 Nonparametric Random Generation
Random observations may be generated from probability distributions by making use of the
fact that the cumulative probability function for any continuous variate is uniformly
distributed over the interval [0, 1] (e.g., Haan, 2002). Thus, for any random variable X with
), the variate FX(X)
29
-------
X
Fx(x)=\fx(u)du (4.8)
-00
is uniformly distributed over [0, 1]. A procedure for generating a random value x fmmfx(x)
starts with generating a random number,/? from the uniform distribution in the interval [0, 1],
then setting R = FX(X), and finally solving for x using the inverse relationship
x = Fxl(R) (4.9)
Unfortunately, an explicit solution to (4.9) is not always possible because FX(X) may not
be a simple function of x. The procedure for solving this equation is achieved in two steps.
First, FX(X) is obtained by substituting the right-hand-side of (4.4) forfx(x) in (4.8) and
commuting the order of summation with the integration operator. The term-wise integration
can be carried either numerically or analytically in terms of resulting error functions.
Secondly, for each randomly generated R, Eq. (4.9) is solved using any of the root searching
techniques. The simple bisection method is used to obtain the root.
A large sequence of independent, identically distributed model residuals wt are randomly
generated using the aforementioned procedure. Latin Hypercube Sampling (LHS) is
employed as an efficient and effective alternative to conventional Monte Carlo sampling
(MCS). The LHS (McKay et al. 1979) divides the CDF (i.e., the F*(x) function) into
segments of equal width from each of which a random variate R is generated. This way the
whole CDF is covered, but with smaller number of replications than by MCS. LHS is more
efficient than MCS at both ends of the CDF. The randomly generated sequence of wt is fed
into Eq. (4.3) to synthesize ensemble time series of st. As indicated at the outset, the
computed st values, in fact, estimate the 3-day average as opposed to daily model prediction
errors.
The generated ensemble of model errors together with the effect of precipitation
uncertainty were used to construct an ensemble of flow duration curves based on a present
land use map (Figure 2a). However, before doing so, model forecast performance was first
examined using the stochastic error model developed above and available streamflow data.
4.3 Model Forecast Evaluation
Figure 7 constitutes the first step in evaluating the model forecast performance during the
post-validation period 5/1/2005-9/30/2005. The simulated streamflows and measured
counterparts are daily and not 3-day averages. The relatively good fit indicated by extending
the validation period further corroborates the calibrated model, and suggests that no further
calibration is warranted. In general, the simulations somewhat overestimated measured
streamflows during high flow periods and underestimated measured values during low flow
periods.
30
-------
To examine model forecast quality, the period from 7/1/2002 to 4/30/2005 (1035 days) is
divided into two parts. The first 900 days (7/1/2002-12/16/2004) are used to construct the
error model and the remaining 135 days from (12/16/2004-4/30/2005) are used for
validation. The synthesized st values are subsequently added to SWAT model simulations
conducted during the same time period to obtain an ensemble of model forecast (500 time
series of streamflows),
ot=Pt+st (4.10)
where the A symbol denotes an estimate. ot is therefore considered as a surrogate to daily
forecast of streamflows. Negative values of ot computed from (4.10) are replaced with zero
flow rates.
Thus, for each t (i.e., day), there are 500 values of computed streamflows,^, available to
construct uncertainty band (confidence interval). Figure 11 compares measured streamflow
with ot. About 7% of observed values fall outside the 95% confidence band. Overall, the
median of simulated 3-day average streamflows compared fairly well with the measured
counterparts. The confidence band tends to be narrower for high flows and wider for low
flows. This indicates that the forecast is more reliable during storm events and least reliable
during smaller events and near base-flow conditions. To further examine model forecast
during low flow conditions, the forecast evaluation period is extended to 9/30/2005. No
significant storm events were recorded during the time period 5/1/2005-9/30/2005. While the
simulated median daily flows generally compare well with measured counterparts, the 95%
confidence band remains relatively wide (Figure 12). It is not clear if this is the result of
inadequate error model fit or poor SWAT model performance during low flow periods. In
general, watershed hydrologic models tend to perform poorer and have larger prediction
uncertainty during low flow events (e.g., Hantush and Kalin, 2005). This phenomenon is
largely attributed to the inadequacy of the models to capture the true intricate nature of
runoff-sub surf ace flow interactions for small events. It will be shown later that simulated low
flows generally have relatively high uncertainty.
-------
60
50
40
0)
D)
| 30
ro
20
10
0
12/16/2004
1/5/2005
1/25/2005
2/14/2005 3/6/2005 3/26/2005 4/15/2005
observed
• (SWAT+error) - median
--•(SWAT+error)-2.5%
•(SWAT+error)-97.5%
Figure 11. SWAT model ensemble forecast (median and 95% confidence band) and
measured streamflows during the validation period (12/16/2004-4/30/2005). A
total of 9 measurements lie outside the 95% confidence band (6.7 %). The
simulated and measured values are 3-day averages.
10
0)
D)
%
&
m
0
5/1/2005
* ' VV ' V i •'' '1JVV V './
5/21/2005 6/10/2005 6/30/2005 7/20/2005 8/9/2005 8/29/2005 9/18/2005
observed
-median
ARMA-normal 2.5%
•ARMA-normal 97.5%
Figure 12. SWAT model forecast (median and 95% confidence band) and measured
streamflows during the post-validation period (5/1/2005-9/30/2005). The
simulated and measured values are 3-day averages.
Since a rationale was set to use the simulated 3-day average streamflows as surrogates to
the daily flows, and for completeness, the ensemble forecast (median and 95% confidence
32
-------
band) and measured daily flows are plotted in Figure 13. The median of 3-day average of the
simulated streamflows appear to have grossly underestimated some of the measured peak
flows. On the other hand, low flows were simulated fairly well.
0 Ll
12/16/2004
1/5/2005
1/25/2005
2/14/2005
3/6/2005 3/26/2005 4/15/2005
observed
-(SWAT+error) - median
•(SWAT+error)-2.5%
•(SWAT+error)-97.5%
Figure 13. SWAT model 3 -day average ensemble forecast (median and 95% confidence
band) and measured daily streamflows during the validation period (12/16/2004-
4/30/2005). A total of 15 measurements fall outside the 95% confidence band (11
In the following subsection we present a methodology for computing SWAT model
response to uncertainty associated with precipitation.
4.4 Uncertainty Due to Precipitation Variability
Precipitation is the central driver of most hydrological processes. Model results are highly
sensitive to precipitation, and uncertainty in precipitation input may affect model output
variability. Other sources of uncertainties, such as the combined effect of parametric
uncertainty, model structural errors and measurement errors, are accounted for by the
stochastic model (4.3). Precipitation temporal variability is addressed here.
For a rational analysis of the potential future impacts of the projected land use alterations
on the hydrologic cycle of the Pocono Creek Watershed, it is essential to generate
meteorological data for future conditions mimicking past conditions. Generated data such as
precipitation, minimum and maximum air temperature, relative humidity, solar radiation and
wind speed need to preserve statistics of past observations. Admittedly, the supposition of
preserving the statistics of the past meteorological data ignores any global or local climate
change effects. Thus, in this analysis, it is assumed that the anticipated changes in the
33
-------
hydrologic regime are solely as a result of the projected modifications in the land use and
land cover. SWAT has a built-in weather generator WXGEN (Sharpley and Williams, 1990)
for atmospheric data generation. WXGEN relies on a first-order Markov-chain model to
define wet and dry periods. For detailed information on weather data generation
methodology, interested readers may refer to Kalin and Hantush (2006b). At the onset of this
study, Mount Pocono station had daily atmospheric data from 10/1/1999 to 6/30/2005, which
means that only 5 to 6 years of data was available to compute monthly statistics. At first
glance, this duration may sound inappropriate for computing climate data statistics as the
SWAT manual recommends the use of 20 years or longer data. Kalin and Hantush (2006b)
investigated the adequacy of relying on a shorter duration data by comparison with rain data
statistics pertinent to the Stroudsburg station, and showed that there are no significant
differences in the statistics of the Stroudsburg station whether the past 5 years or 20 years of
data is used. Therefore, for this specific application, it may be concluded that 5 years
precipitation data reasonably reflects the statistics of a 20 years record for Pocono Creek.
4.5 Monte Carlo Simulations
To take into account the precipitation uncertainty, WXGEN is used to generate 500 sets of
20-year long records of daily precipitation assumed to represent precipitation from 1/1/2005
to 12/31/2024 based on averaged historical rainfall statistics at the Mount Pocono and
Stroudsburg stations. A warm-up period from 1/1/1975 to 12/31/2004 of 30 years is used to
eliminate the effect of uncertain initial conditions. Each of the synthesized time series of
daily precipitation is fed into SWAT with the current land use map (Figure 2(a)) to obtain a
20-year long time series of daily streamflows. For each MC simulation, 3-day moving
average is computed from the simulated daily streamflows. The MC simulation is repeated
for a total number of 500, 20-year long time series of daily streamflows, from which the
ensemble of 3-day averaged flows are computed.
Assuming that future simulations errors preserve the statistical characteristics of the
historical record and follow the same stochastic model obtained in Eq. (4.3), an ensemble
forecast of streamflows (median and 95% confidence band) for the period 1/1/2005 to
12/31/2024 is constructed by generating an ensemble of 20-year long sequences of st which
then are added to the ensemble of MC simulations. In this manner, all sources of errors are
accounted for in the forecast. The sequence of independent, identically distributed wt are
randomly generated according to the procedure outlined in section 4.2.3, then used in
conjunction with Eq. (4.7) to synthesize an ensemble of 500 20-year long sequences of daily
(actually 3-day average) st values. The MC ensemble of 20-year long sequences of daily
streamflows generated above are then corrupted by the randomly generated errors (st) to
obtain the ensemble forecast of streamflows at the USGS gauge station according to Eq.
(4.10). It is worthwhile to note that the streamflows computed according to (4.10) constitute
what is expected to be observed or measured rather than actual flow rates.
Summary statistics of some of the model outcomes for the current land use scenario are
presented in Table 5. Mean, standard deviation (std), coefficient of variation (CV), median,
and 2.5th and 97.5th percentiles (95% confidence interval) of the 500 realizations for each
34
-------
flow statistics are given in the table to explicate uncertainty involved in the model outcomes
due to precipitation as well as structural, parametric, and observation errors.
Table 5. Summary Statistics of Ensemble Streamflow Forecast Based on Most
Recent land use (Figure 2a). Values in the Table are Obtained From
500 Time Series of Synthesized Daily Streamflows
Average Daily Average Monthly Average Monthly Average Annual
Flow (m3/s) Median Daily Maximum Daily Maximum Daily
Flow (mVs) Flow (m3/s) Flow (mVs)
(1) (2) (3J (4)
mean 2.42 2.10 7.12 17.39
std 0.203 0.182 0.411 1.93
CV 0.084 0.087 0.058 0.111
median 2.42 2.10 7.13 17.43
95% C.I [2.02,2.80] [1.75,2.45] [6.35,7.91] [13.91,21.82]
NOTE: Daily flow here refers to 3 day moving average flow.
In the table, summary statistics related to high flow conditions as well as long term
averages are given. The term "Average" in each column title denotes arithmetic average over
the 20-year simulation period. For example, for each realization (i.e., time series) out of 500,
daily flows are averaged over the 20-year simulation record to yield "average daily flow".
Thus, there are 500 such "average daily flow" values from which the mean, median, CV, and
95% confidence limits are computed. Average monthly-median of simulated daily
Streamflow in column (2) represents the arithmetic average over the simulation period of the
median daily stream flows for each month and as indicated above is a required input to the
PIFM wild brown trout habitat model. In column (3) the monthly maximum daily
streamflows averaged over the simulation period are given. Column (4) shows the annual
maximum daily streamflows averaged over the simulation period. Note that the summary
statistics in Table 5 assume that the watershed remains undisturbed in the next 20 years,
which obviously is contrary to what is anticipated. As such, they are hypothetical and are
intended solely to provide insights into the effect of uncertainties on the interpretation of
model predictions. Nevertheless, assuming that the land use/landscape in the watershed did
not undergo significant changes in the past twenty years, the tabulated results may provide
insights into the current flow conditions in Pocono Creek Watershed. The following remarks
can be made from Table 5:
35
-------
i. The average daily flow, average monthly median of daily flows, and average monthly
maximum daily flow have small CV values and relatively narrower 95% confidence
bands. Given all sources of uncertainty, these flow measures show small variability and
appear to be reliable measures of flow characteristics.
ii. The average annual maximum daily flow shows a relatively greater uncertainty, as the
larger CV and wider 95% confidence band indicate. The performance of flood control
measures designed based on annual maximum daily flows should factor in the
computed uncertainty.
iii. The low uncertainty associated with median monthly daily flow is rather encouraging
given that it is an important parameter for the wild brown trout habitat (PIFM) model.
The ensemble characteristics of the simulated flow duration curves are depicted in Figure
14. The figure plots the median of the flow duration curves and the lower (2.5th percentile)
and upper (97.5th percentile) limits of the 95% confidence interval. Each of the 500
synthesized o^time series yielded one duration curve. A flow duration curve is a plot of the
flow rate (say, x) versus probability of exceedance, P(X> x). The return period, T, defined as
the average recurrence interval between events equaling or exceeding a specified flow
magnitude, x, is the reciprocal of P(X > x): T = II P(X> x). High flows are associated with
small probability of exceedance, whereas low flows are associated with large probability of
exceedance. The ensemble of flow duration curves in Figure 14 can be interpreted as follows.
From the figure it can bee seen that for P(X > x) = 0.002, daily streamflow at the gauge
station is between 16 and 25 mVs with 95% confidence. Since the corresponding return
period is T = 1/0.002 = 500 days, this means that 500 day return period flow is between 16
and 25 m3/s with 95% confidence.
One may be lured into believing that uncertainty increases with increasing daily
streamflow as the 95% confidence band gets wider with decreasing probability of
exceedance, however, careful inspection of Figure 15 reveals that low flows have greater CV
values and, thus, show much higher uncertainty than medium range and high flows. It is
interesting to note that predicted streamflow rates ranging from about 2 to 11 (m3/s) daily
flow rates, corresponding to return periods of 2.5 and 50 days, respectively, have the smallest
CV values and, therefore, lowest uncertainty. Higher flows with a recurrence period longer
than 50 days have relatively higher uncertainty but much smaller than that associated with
low flows, with return periods smaller than 2 days. Of course, these results are specific to
Pocono Creek watershed.
36
-------
median
95% C.I- upper limit
95% C.I-lower limit
0.001
0.01 0.1
Probability of exceedance
Figure 14. Median and 95% confidence band for the MC simulated duration curve. The flow
duration curves are generated from 500 replicas of 20 years daily streamflows,
c n A
0 U'4
'ra
-2 03
"o
"c
m n o
Coeffici
o c
-^ K
n
O.C
^,
01
"^
^
^
-
t~
0
'"^~~1
01
—
-
I
~\~\~
!
0.1
—
-
r
~[~i
^
1
Probability of exceedance
0.01 0.1
1-(Probability of exceedance)
Figure 15. Coefficient of variation computed from ensemble of flow duration curves versus
probability of exceedance in the left panel P(X > x), and cumulative probability
P(X< x) in the right panel.
4.6 Conclusions
In this section model prediction capability was investigated using time series analysis and
Monte Carlo-type Latin Hypercube simulations. Structural errors, parametric uncertainty,
measurement errors, and rainfall variability together contribute to the SWAT model
37
-------
predictive uncertainty. To minimize the effect of systematic SWAT prediction errors during
significant events, 3-day averaged streamflows were considered as surrogates to daily flows.
A stochastic model, which lumps the first three sources of errors, was developed by fitting
ARIMA (1,1,1) x (0,0,1)3 time series model to 3-day averaged observed daily model errors
during the calibration and verification period. A split-sample approach was implemented to
construct and validate the error model. The stochastic error model allowed for constructing a
forecast band (the median and 95% confidence band) for the SWAT model simulations. It
was shown that most of the observed streamflows during the validation period (12/16/2004-
4/30/2005) were within the 95% confidence band. The constructed model forecast was
further evaluated using an additional set of measured streamflows during the period
5/1/2005-9/30/2005. Although the results showed consistency, the 95% confidence band was
relatively wide due to the relatively low-flow period.
The effect of temporal precipitation variability was also accounted for using historical
precipitation data recorded at the Mount Pocono and Stroudsburg stations. Synthesized
precipitation rates by the SWAT built-in weather generator WXGEN along with MC
simulation of the SWAT model and generated model errors all together were used to
construct an ensemble forecast of daily streamflows for the next twenty years. It was shown
that given all sources of model uncertainty, the average daily flow, average monthly median
of daily flows, and average monthly maximum daily flow could be reliably predicted. The
relatively low uncertainty associated with monthly median of daily flows indicate that this
flow measure may be used reliably as an input to the brown trout habitat model (PIFM).
Averaged over the simulated 20-years period, annual maximum daily flow showed relatively
greater variability. The ensemble of daily flow duration curves was summarized by the
median and 95% confidence band. The ensemble duration curve graph allows for the
estimation of the daily flow range with 95% confidence for a specified design recurrence
(return) period. SWAT simulated daily streamflow rates ranging from about 2 to 11 (m3/s)
showed least uncertainty. Computed daily streamflow rates below 2 m3/s had the greatest
uncertainty, whereas for higher than 11 m3/s, uncertainty was moderate. The higher
uncertainty associated with low flows was not surprising as watershed models tend to
perform relatively poorer during small events.
38
-------
5 Impact of Land Use Changes
5.1 Introduction
Over the past 25 years, the population of the United States has grown over 30% (USDC
Census Bureau, 2005). Naturally, such an ample growth in population leads to substantial
increase in urbanized areas and results in a degradation and loss of forested and agricultural
lands. Hydrologically, urbanization is accompanied by increased imperviousness in the
landscape. Increase in impervious areas as well as reduction in soil permeability results in
reduced infiltration rate, which means that more precipitation becomes surface runoff and
less water is recharged to ground water. When increased surface runoff is combined with the
effect of reduced surface roughness - a consequence of which is a shorter travel time - it is
inevitable to observe more frequent and more intense local flood events. Such alterations in
the flow regimes of streams and channels may lead to changes in channel morphology in the
form of channel widening and deepening. On the other hand, groundwater recharge reduction
results in the drop of groundwater levels and reduction of base flow to stream flow.
Groundwater depletion can be further amplified due to an increase in groundwater
withdrawal that accompanies population and economic growth. Consequently, stream flow is
further depleted as groundwater levels are increasingly lowered by increased pumping. This
could have undesirable ecological consequences not only due to limited available water, but
also because of increase in pollutant concentrations and limited capabilities of the streams to
dilute any toxic spills.
As indicated earlier, Monroe County, where the Pocono Creek is located, has the second
fastest growing population in the state of Pennsylvania. By the year 2020, a
60% increase in Monroe County population is forecasted, and it is further projected that more
than 70% of the Pocono Creek watershed will turn into commercial and residential areas,
currently standing around 6% (Figure 16). The impact of population growth and urbanization
on potential alterations in surface and groundwater regimes will be assessed through the
modeling framework that has been developed.
5.2 Present Land Use Map and Future Build Out
A simulation period of 20 years, from 1/1/2005 to 12/31/2024, is employed to study the
impact of changes in the land use on the hydrology of the watershed. Two land use scenarios
are considered. First scenario (LU2000) assumes that land use pattern of the year 2000, as
depicted in the top panel of Figure 16, is preserved during the simulation period. In the
second scenario (LU2020) it is assumed that the build out would occur in 2020 or after
(bottom panel in Figure 16). In both scenarios, land use pattern is assumed to remain the
same throughout the 20 year simulation period. In other words, land use pattern is assumed
time invariant during the course of the simulations. Ideally, it is desirable to let the land use
pattern change gradually during the course of the simulations. However, this is not an easy
task due to model limitations as this will require dynamic updating of land use related model
parameters with time. Recall that the watershed model was calibrated and validated based on
the land use pattern LU2000. The model setup for simulation of LU2000 is extended to
39
-------
perform the model simulations of LU2020.
watershed
Open Water
Perennial Ice/Snow
Low Intensity Residential
High Intensity Residential
Commercial/Industrial
Bare Rock/Sand/Clay
Quarries/Strip Mines
Transitional
Deciduous Forest
Evergreen Forest
Mixed Forest
Shrubland
Orchards/Vineyards
Grasslands/Herbaceous
Pasture/Hay
Row Crops
Small Grains
Fallow
Urban/Recreational Grasses
Woody Wetlands
Emergent Herbaceous Wetlands
LAHD USE
Residential-High Density
Residential-Medium Density
Residential-Low Density
Conmer c ial/Transportat ion
Pasture
Water
Wetlands
Forest
Agricultural Land-Row Crops
TOTAL
2000
0.05
3.53
2 .25
3 .52
43
80
1
3
85.23
0. 19
2020
0.77
8.00
22 .84
0.27
1.43
3.80
18.70
100.00 100.00
Figure 16. Distribution of land use pattern in the Pocono Creek watershed in year 2000 (top)
and the projected land use pattern for the year 2020 (bottom).
In simulating the effect of the two land use scenarios on streamflow characteristics all
sources of model output error are neglected except rainfall variability; i.e., only SWAT
model output bt is considered, with the ensemble of rainfall rates generated for the period
from 1/1/2005 to 12/31/2024. This is done primarily for two reasons. First, the error
stochastic model (4.3) applies to LU2000, and it is therefore not conclusive if a similar error
structure would apply to model simulations based on LU2020. Secondly, it is the relative
40
-------
changes in the streamflow characteristics due to projected land use changes that are of
concern rather than absolute predictions. Computed relative changes tend to have
significantly lower uncertainty, consequently, they are much more reliable than absolute
predictions of daily streamflows.
5.3 Monte Carlo Simulation
To take into account the precipitation uncertainly, we generated 50 sets of distinct daily
precipitation records of 20 years length, each assumed to represent precipitation from
1/1/2005 to 12/31/2024. Measured precipitation data from 1/1/1975 to 12/31/2004 is inserted
at the beginning of each record to obtain 50 precipitation input data files, each of which
contain 50 years of daily precipitation. For each scenario, model simulations are performed
for each of the 50 precipitation data files. A total of 2x50x50=5000 years of model
simulations are performed at the daily time scale. Again, the first 30 years (1/1/1975 to
12/31/2004) of each realization is ignored for model warm-up purposes. Only the last 20
years of each realization are retained for further analysis. The reason for using the same
measured precipitation data during the 30 year warm-up periods of all realizations is to
minimize uncertainty relevant to initial conditions. Experimentation revealed that 50
realizations were adequate, and that more realizations had a negligible effect on the model
outcome.
5.4 Predicted Changes in Streamflow
The MC simulations yielded 50 time series for each of the SWAT Model outputs that
include, among others, stream flow, base flow, and groundwater recharge. Summary statistics
of some of the model outcomes for the two scenarios, LU2000 and LU2020, are presented in
Table 6 along with relative changes that might be expected when the land use pattern in the
Pocono Creek Watershed changes from the one given in year 2000 to the one projected past
year 2020. The statistics in the far left column and associated results in the table correspond
to the variations of model output. Mean, standard deviation (std), coefficient of variation
(CV), median, and 5th and 95th percentiles (90% confidence interval) for each design flow are
computed from the resulting MC simulation and are given in the table to explicate
uncertainty involved in the model outcomes due to precipitation. In the table summary
statistics related to low and high flow conditions as well as long term averages are given. The
first two columns given in the table are essentially base flow (BF) and stream flow (SF)
averaged over the 20-year simulation period. The lowest computed flow occurring once
every 10 years averaged over a 7-consecutive-day period (7Q10) listed in column (3) is
widely used as a low-flow index in the United States (Smakhtin, 2001). This hydrologically-
based design flow parameter is also used to protect against chronic effects by requiring that
water quality criteria must be met at all times except during the 7Q10. Average monthly-
median of simulated daily SF in column (4) represents the arithmetic average over the
simulation period of the median daily stream flows for each month and as indicated above is
a required input to the PIFM wild brown trout habitat model. In column (5) and (6) the
monthly and annual maximum daily stream flows averaged over the simulation period are
given, respectively.
41
-------
Table 6. Summary Statistics of Computed Streamflow Characteristics at the USGS Gauge
Station. The Results are Derived From 50 MC Simulations Each 20 Years Long
mean
std
CV
median
[5, 95]%
LU2000
LU2020
% change
LU2000
LU2020
LU2000
LU2020
LU2000
LU2020
% change
LU2000
LU2020
% change
Average BF
(m3/s)
(1)
1.314
0.911
-30.6%
0.071
0.056
0.054
0.062
1.318
0.909
-31.0%
[1.205,1.421]
[0.828,0.993]
Average SF
(m3/s)
(2)
2.452
2.426
-1.1%
0.137
0.140
0.056
0.058
2.457
2.431
-1.1%
[2.232,2.682]
[2.200,2.664]
7Q10
(m3/s)
(3)
0.284
0.252
-11.1%
0.054
0.055
0.191
0.218
0.282
0.242
-14.4%
[0.197,0.364]
[0.167,0.326]
Average
Monthly
Median
Daily SF
(m3/s)
(4)
2.033
1.829
-10.1%
0.108
0.104
0.053
0.057
2.027
1.826
-9.9%
[1.852,2.207]
[1.650,2.002]
Average
Monthly
Maximum
Daily SF
(m3/s)
(5)
6.11
7.40
21.1%
0.450
0.526
0.074
0.071
6.17
7.47
21.0%
[5.35,6.71]
[6.52,8.17]
[21.8,21.7]%
Average
Annual
Maximum
Daily Flow
(m3/s)
(6)
17.47
20.86
19.4%
2.28
2.67
0.131
0.128
17.40
20.81
19.6%
[13.59,22.31]
[16.26,26.55]
[19.6,19.0]%
mean: Arithmetic average of 50 realizations
std: Standard deviation of 50 realizations
CV: Coefficient of variation (std/mean)
median: 50th percentile of 50 realizations
[5,95]%: 5th and 95th percentiles of the 50 realizations
BF: Base flow contribution to stream flow
SF: Stream flow
7Q10: Seven days average low flow with a 10 year
return period.
LU2000: Simulations performed using land use
coverage for year 2000
LU2020: Simulations performed using land use
projections for 2020
% change: (LU2020-LU2000)/LU2000
Table 1 reveals the following observations and associated conclusions:
i. The mean value of 50 realizations for the average of simulated BF is expected to decline
by 31%. This is simply due to reduced infiltration and consequently less recharge to the
ground water. This reduction can be further aggravated if the effect of increased
groundwater withdrawal due to higher water demand associated with population growth
is imposed. To explore the effect of potential increase for groundwater demand, a more
detailed groundwater flow model (MODFLOW) is being developed for the region by the
USGS. In the simulations, groundwater pumpage is ignored. Precipitation uncertainty
does not seem to impart much uncertainty on BF as CV is less than 0.1 for both
42
-------
scenarios. The median, 5th and 95th percentiles are also expected to decline by the same
rate as the mean.
ii. The average of simulated daily SF is relatively unaffected by urbanization. This can be
explained on the basis that reduction in BF is offset by the increase in surface runoff. Yet,
a negligible reduction (1.1%) in SF is observed. Mass balance of soil water dictates that
average evapotranspiration (ET) is not impacted much by land use changes. Indeed, close
examination of ET outputs from the model confirms this. Although forested areas have
high transpiration (T), impervious surfaces have higher solar reflectivity, and thus higher
evaporation (E). The reduction in (T) seems to be compensated by the increased (E).
Similar to BF, uncertainty in precipitation has insignificant impact on SF. It should be
noted that since the average of simulated daily SF is obtained by averaging simulated
daily streamflows over the entire simulation record, it therefore does not represent actual
daily flow conditions.
iii. On average, the computed 7Q10 is expected to decrease by almost 11%. This reduction
can be attributed to the reduction in BF as low flows in perennial streams, such as Pocono
Creek, are driven by BF which is predicted to decline by 31% on average over the study
watershed. With 90% confidence, 7Q10 can be expected to be within the interval [0.197,
0.364] m3/s with the current land use conditions. With the projected land use conditions
the 90% confidence interval for 7Q10 is [0.167, 0.326] m3/s. The CV is estimated to
increase by 14% due to changes in land use leading us to conclude that the combined
effect of precipitation uncertainty and projected land use alterations are likely to cause
higher uncertainty in 7Q10.
iv. Mean, median, 5th and 95th percentiles of average (i.e., averaged over the simulation
period) monthly median of simulated daily SF as shown in column (4) are all reduced by
about 10%. It is also clear from the table that the levels of uncertainties in both scenarios
do not differ considerably, as they are only marginal. From a hydrologist's perspective, it
is hard to reach a conclusion on how much risk this reduction poses to the brown trout
population in Pocono Creek. Simulations based on the PIFM habitat model and the
computed median flows will provide answers.
v. As expected, average monthly maximum of simulated daily SF increases as summarized
in column (5). The expected increase is around 21% for mean, median, and 5th and 95th
percentiles of the 50 realizations. Reasons for increase in peak flow magnitude as well as
peak flow frequencies as a result of urbanization are well known: reduced infiltration rate
and shorter travel times. The levels of uncertainties are about the same for both scenarios
and can be deemed negligible with CV values of 0.07.
vi. Similarly, the increase in the average of annual maximum daily flow is anticipated with
LU2020 projections. The mean and the median of average of annual maximum daily flow
are predicted to increase by about 19%. The uncertainty associated with this design flow
is the second highest after 7Q10, with a CV greater than 0.1. The 90% confidence band is
estimated to increase from [13.59, 22.31] to [16.26, 26.55] in units of mVs.
In Figure 17, spatial distributions of the 20-year simulation average of annual
groundwater recharge estimates [mm] at the subbasin scale are shown for the two scenarios.
These are only model estimates, as real measurements of groundwater recharge are not
43
-------
available. The figure delivers some important insights. There is a significant spatial variation
of groundwater recharge in the watershed. Lower portions of the watershed (associated with
lower elevations) have higher recharge rates than upper portions. When recharge estimates of
the scenarios are compared, it is seen that with few exceptions almost all sub watersheds
experience reduction of groundwater recharge. Another important observation is related to
the extent of spatial variation of recharge in the two scenarios. In LU2000, the computed
average recharge rates vary from 155 mm/year to 576 mm/year. In LU2020 the computed
range of groundwater recharge is from 131 mm/year to 432 mm/year. This indicates that the
projected land use change is expected to slightly reduce the spatial variation of groundwater
recharge over the Pocono Creek watershed.
LU2000
Annual Recharge (mm)
| [100-150
| [150-200
| 200 - 250
j 250 - 300
| 300 - 350
] 350 - 400
] 400 - 450
I 450 - 500
]500 - 550
I 550 - 600
Figure 17. Annual groundwater recharge distributions in the Pocono Creek watershed for the
two land use scenarios.
It should be noted that the results in Table 6 and Figure 17 assume that the build out in
LU2020 occurs by 2020. The results would still characterize the watershed response even if
the build out in LU2020 is presumed to occur after 2020. This may be true provided that
pattern and variability of future precipitation remains unaltered.
The predicted annual groundwater recharge distributions depicted in Figure 17 were
tabulated and passed on to USGS (Malvern, PA) to calibrate a quasi-three dimensional
groundwater flow model (MODFLOW) for the watershed. The computed monthly median of
daily flows were tabulated and passed on to the PA F&B Commission for use as input data to
the brown trout habitat model (PIFM). Annual groundwater recharge averaged over the
simulation record and spatially over the watershed was predicted to decline by 31%, from
424 to 292 mm/yr based on the LU2020 projections. The coefficient of variation (CV) of the
computed annual groundwater recharge did not exceed 0.08, indicating fairly reliable model
estimates thereof.
44
-------
5.5 Predicted Changes to Flow Frequency and Duration
The medians of daily duration curves based on LU2000 and LU2020 are plotted in Figure 18.
The probability of exceeding high flows and the risk for flood hazard are predicted to
increase based on LU2020. On the other hand, flow exceeded at least 90% of the time (< 0.8
m3/s) decreases in moving from LU2000 to LU2020, as depicted by the inner panel (Figure
18). Equivalently, this means that the chance for the flow to be less than or equal to a given
low flow threshold increases. Therefore, base flows are likely to decrease with the projected
land use changes. This is, of course, excluding the effect of the anticipated increase in
groundwater withdrawals which will further decline base flows. Also note from Figure 18
that median flow due to LU2020 is slightly smaller than the one due to LU2000.
30
25
I/T 20
r)
E
I 15
90% X95% 100%
001
0.01 0.1
Probability of exceedance
Figure 18. Median of the ensemble of the MC-simulated duration curves, bold thick line
corresponds to simulations with LU2020, thin line corresponds to simulations
with LU2000. The flow duration curves are generated from 50 replicas of 20
years of daily streamflows, ot.
5.6 Conclusions
The Pocono Creek watershed is threatened by high population growth anticipated over the
next two to four decades. Potential effects of population growth and urbanization have
heightened the need for implementing sustainable water resource management strategies in
the watershed. In this section, potential hydrological changes in Pocono Creek due to
anticipated build out in the watershed were investigated using the calibrated SWAT model
and Monte Carlo simulations. Simulations accounted for anticipated rainfall variability over
a 20 year period, with the current land use pattern and projected build out for the year past
2020, named LU2000 and LU2020, respectively. Simulation results revealed that on the
average, daily base flow is expected to be reduced by 31%. However, this is not expected to
45
-------
cause a noteworthy reduction in average daily streamflow (averaged over the entire
simulation period) as SWAT predicted that reduction in base flow would be balanced by
increased surface runoff. The computed low-flow index, 7Q10, is expected to decline by 11%
due to anticipated base-flow reduction. A metric for the sustainability of fish habitat, the
computed monthly median daily flow is expected to decline by 10% on the average. The
monthly peak of simulated daily flows and annual maximum daily flow are expected to
increase by 21% and 19%, respectively. The computed 7Q10 and average annual maximum
daily flow showed relatively higher uncertainty than the other flow characteristics. Projected
build out in the watershed is estimated to cause a significant decline in the annual
groundwater recharge rates. The decline in the watershed-averaged annual groundwater
recharge is predicted to be 31% based on LU2020 projections. The spatial variability of
groundwater recharge appears to diminish with the projected urbanization in the watershed.
In general, the likelihood that the watershed will experience high and low streamflows will
increase with the projected urbanization, as indicated by the median of the MC simulated
flow duration curves.
46
-------
6 Critical Source Areas
6.1 Introduction
It is now evident from model predictions that the projected land use changes in the Pocono
Creek watershed have the potential of increasing average annual maximum and average
maximum monthly flows, and reducing 7Q10 and average median monthly daily flows at the
watershed outlet. Informed management decisions may benefit from the identification of
portions of the watershed that have the highest contribution to the reduction/increase in the
quantity of interest. In a sense, preserving the land use of a particular area in the watershed
can be considered as a best management practice (BMP). From this point of view, the
problem can also be posed as identifying the locations of those BMPs that minimize the
predicted changes. Of course, socioeconomic and policy matters may interfere with and
preclude the implementation of a particular BMP, but knowing beforehand (i.e., before land
development) critical areas in the watershed provides science-based guidance to the planning
process.
There are limited applications in literature looking at the aspect of identifying or
apportioning areas within a watershed based on their relative contributions to flow at the
outlet. Saghafian and Khosroshahi (2005) address the same problem by focusing on flood
source areas and they also emphasize lack of applications of this type. The approach adopted
here, to be addressed shortly, has some similarities, however it differs from their work in two
ways: i) in addition to high flows, the focus is extended to low and median flow
characteristics (7Q10, monthly median daily flows), and ii) the focus here is not on
individual event hydrographs, but rather on continuous flow time series.
6.2 Methodology
Let us assume that a given watershed is divided into k number of subareas, which will be
referred as "elements" henceforth. There are two land use conditions, current (c) and future
(f). The model run at the daily time scale with the future land use scenario for a given
duration generates the flow time series Qf at the watershed outlet. Suppose that all the
elements have the future land use/cover with the exception of element y retaining its present
status. The generated flow time series at the outlet with this setup will be denoted as Qf'J .
For both flow time series, Qf and Qf'J , the flow characteristics of interest, say F is
computed and designated as Ff and FfJ respectively. The following two indexes are
defined to assess the relative impact of element y' on the flow characteristic F:
i i
Ff AjlA*
where Aj and Aw indicate areas of element y' and the whole watershed, respectively. The first
index, a/, signifies the absolute impact of element y on F. The second index, /^, which is
basically normalization of q with the percentage area of element y, measures the impact of
land use changes in the element on F assuming these very changes occur over the entire
47
-------
watershed. In other words, /?is suitable for assessing the impact of land use changes per unit
area. By computing a/ and $ for j = l,..,k one can rank the areas in the watershed from most
critical to least. It is important to note that the order of ranking could be different not only
depending on the index used but also depending on the flow characteristic of interest, F.
6.3 Results
As part of the modeling requirement, the Pocono Creek watershed was divided into 29
catchments, numbered from 1 to 29 with catchment 29 being the most downstream (Figure
19a). For management purposes this discretization might be too detailed. Hence, as an
alternative we divided the watershed into 7 larger areas which are combinations of the 29
subbasins (Figure 19b). It should be noted that the larger areas correspond to the
management areas in the pilot study (2001), with areas 4,5, and 6 in Figure 19(b) together
falling within the management area 3 of the pilot study (Figure 20). Area 7 in Figure 19(b)
overlaps with management areas 2 and 3 of the pilot study (Figure 20).
(a)
Figure 19. Watershed subdivisions used in determination of critical areas, (a) finer, (b)
coarser.
48
-------
Figure 20. Management areas in Pocono Creek Watershed according to Pocono Pilot Study.
Table 7 summarizes the computed indices, a and (3 with the coarser spatial scale for
7Q10, monthly median of daily flows and annual maximum daily flows, which divides the
watershed into 7 drainage areas, j=l,..,7. Also given in the table are the rankings of each area
for each index and flow characteristics. Figure 21 depicts these rankings on a gray scale
spectrum for visual purposes. One immediate observation is that the rankings based on a and
/? are quite consistent for monthly median of daily flows and annual maximum daily flow.
Areas having the highest and lowest impact on these two flow characteristics do not change
with the type of index. Management area 4 has the highest impact on the reduction of 7Q10
when a index is used, whereas management area 5 is the biggest contributor when /? index
issued.
Table 7. Computed Indices a and /? for the 7 Management Areas for 7Q10, Monthly
Median of Daily Flows and Annual Maximum Daily Flows
rank
1
2
3
4
5
6
7
7Q10
j
4
7
2
5
1
6
3
a\
8.58%
7.15%
6.18%
4.95%
3.60%
1 .95%
1 .47%
I
5
4
7
2
6
3
1
i.
0.54
0.51
0.41
0.35
0.22
0.17
0.17
monthly median of daily flow
I
4
7
2
1
5
6
3
a\
2.55%
2.51 %
2.03%
1 .97%
0.85%
0.73%
0.59%
I
4
7
2
5
1
6
3
&
0.15
0.14
0.12
0.09
0.09
0.08
0.07
annual maximum daily flow
i
7
4
1
6
2
3
5
a\
-5.48%
-3.88%
-2.47%
-1 .96%
-1 .69%
-0.56%
-0.50%
I
7
4
6
1
2
3
5
i
-0.32
-0.23
-0.22
-0.11
-0.10
-0.07
-0.05
49
-------
Figure 21. Ranking of the 7 management areas based on a (a, b, c) and /?(d, e, f) indices for
7Q10, monthly median of daily flows, and annual maximum daily flow,
respectively.
Based on a, area 4 (Figure 19(b)) is ranked first in terms of impact on 7Q10 and monthly
median of daily flows, whereas areas 6 and 3 ranked the lowest, 6th and 7th, respectively.
Comparison of Figure 21 with Figure 17 reveals that area 4, through which Bisbing Run and
Bulgers Run flow (Figure 20), contributes the highest annual groundwater recharge among
the seven areas depicted in Figure 19(a). Since base flow is directly related to groundwater
recharge, it is therefore expected that area 4 would have the highest impact on 7Q10. The
predicted significant reduction in the annual groundwater recharge in area 4 due to projected
land use changes (Figure 17) means a significantly greater fraction of precipitation would be
available for runoff. This may explain area 4 being ranked first in terms of impact of land use
changes on monthly median of daily flows, and ranked second for annual maximum daily
flows.
Area 7 is ranked second based on a for its impact on 7Q10 and monthly median of daily
flows, and first for annual maximum daily flows. Close proximity to the USGS streamflow
gauge station (i.e., being downstream most), relatively high groundwater recharge rates, and
the relatively larger area may be contributing to this ranking.
Significant portion of area 6 is occupied by wetlands (about 50%), which are maintained
in both LU2000 and LU2020 as flow through features with zero infiltration. It is assumed
that wetlands are preserved during land development. Therefore, groundwater recharge in
area 6 is slightly reduced, thus, leading to least impact on 7Q10 and lowest ranking based on
a. Area 3, where the Camel Back Ski Area is located and Coolmoor Creek flows (Figure 20),
is characterized by steep topography that is conducive to high runoff and low infiltration
50
-------
regardless of the landscape features. This explains the lowest and second lowest ranking of
area 3 with a.
The ranking based on /? index is useful in assessing impacts when only portions of each
of the 7 areas are planned to be developed. The most notable is area 5 which is ranked first
for 7Q10 using /? index. With the exception of a small, steep northern portion, the area is
topographically characterized by relatively low steepness and long slope length (Figure 1),
and contributes significant groundwater recharge (Figure 17). Developments planned
according to LU2020 may have the highest impact on 7Q10 if they occur within area 5.
6.4 Conclusions
The Pocono Creek watershed was divided into seven catchment areas and an index
methodology was presented to rank the catchments' areas based on their impact on key
streamflow characteristics due to anticipated land disturbances. The catchments' areas
conformed to the six management areas in the pilot study (2001). Two indices were
proposed. The first index, a, signifies the absolute impact of a particular catchment area on
the watershed response. The second index, /?, is a normalized by the percentage area of the
catchment, and therefore describes the impact of land use changes per unit area.
Using these indices, catchments' areas were ranked based on their potential to impact
changes in the streamflow characteristics due to projected build out in the watershed. With a
few exceptions, a and |3 indices produced similar rankings among the 7 catchment areas for
7Q10, monthly median of daily flow, and annual maximum daily flow. Groundwater
recharge, area, topographic features, and proximity to the streamflow gauge station may have
contributed to the ranking results. The most downstream catchment, area 7, ranked first in
terms of impact on annual maximum daily flows, and second in terms of impact on 7Q10 and
monthly median daily flows. Catchment area 4, associated with the highest groundwater
recharge, was ranked first and second based on a and |3 indices, respectively, with regard to
impact on 7Q10. In general, areas characterized by steep topography and significant wetlands
ranked low, some times the lowest, with respect to impact on changes in the three design
flows.
51
-------
7 Summary and Conclusions
The Pocono Creek watershed, which is located in Monroe County, PA, is threatened by high
population growth and urbanization. Of concern specifically is the potential impact of future
developments in the watershed on the reduction of base flow and the consequent risk for
degradation of wild brown trout habitats in Pocono Creek. Anticipated increase in
imperviousness, on the other hand, is expected to elevate the risk for floods and the
associated environmental damage. A watershed hydrologic modeling study was initiated by
the U.S. EPA in collaboration with the U.S. Geological Survey and the Pennsylvania Fish
and Boat Commission to assist Monroe County in planning for sustainable future
developments in the Pocono Creek watershed.
Good application track-record of the SWAT model and its suitability to address
watershed management problems led to the selection of the model to achieve the objective of
quantifying the impact of anticipated land use changes on the hydrologic response of the
Pocono Creek watershed. The model was successfully calibrated and validated for two
sources of precipitation data, raingauge measured and radar estimated (NEXRAD). Two
raingauge stations were available outside the watershed, but in close proximity to the
perimeter. Simulated daily and monthly streamflows at the outlet compared fairly well to the
observed values, for both sources of the rainfall data. The results reinforced the notion that
NEXRAD is an effective source of spatio-temporal precipitation data and a viable alternative
to the very costly installation, operation, and maintenance of a raingauge network. Future
modeling studies in ungauged watersheds can be conducted with NEXRAD rainfall data.
Stochastic error propagation analysis was conducted using time series analysis and MC
simulations to evaluate model predictive uncertainty. Model errors accounted for model
structure uncertainty, parametric uncertainty, measurement errors, and rainfall variability. It
was shown that the calibrated model was consistent in its forecast capability. MC simulations
over a 20-year long period yielded an ensemble of rating curves of which the median and
95% confidence band of daily streamflows were plotted. These plots allow for the
construction of the 95% confidence band for design flows corresponding to any given
recurrence or return period. SWAT simulated daily streamflow rates in the range 2 to 11
(m3/s) showed least uncertainty. Computed daily streamflow rates below 2 m3/s had the
greatest uncertainty, whereas for higher than 11 m3/s uncertainty was moderate.
MC simulation over a 20-year period showed that the average daily base flow is expected
to be reduced by 31% should the projected build out in the watershed occur. However, this is
not believed to cause a noteworthy reduction in average daily streamflow as the model
predicted that reduction in base flow would be balanced by increased surface runoff. The
computed low flow index, 7Q10 is expected to decline by 11%, and the monthly median
daily flow is expected to be reduced by 10% on the average. Monthly peak of simulated daily
flows and annual maximum daily flow were predicted to increase by 21% and 19% on the
average, respectively. Groundwater recharge rate averaged over the watershed was predicted
to decline by 31% due to the projected land use changes. The median of the MC simulated
flow duration curves showed that, in general, the likelihood that the watershed will
experience high and low streamflows will increase with the projected urbanization.
52
-------
An index methodology was developed to rank seven subwatersheds - composing the
modeled area of the Pocono Creek watershed - based on their relative impact on the
watershed response to land developments. The first index, a, signifies the absolute impact of
a particular catchment area on the watershed response. The second index, /?, is a normalized
by the percentage area of the catchment, and therefore describes the impact of land use
changes per unit area. With a few exceptions, a and |3 indices produced similar rankings
among the 7 catchment areas for 7Q10, monthly median of daily flow, and annual maximum
daily flow. Groundwater recharge, area, topographic features, and proximity to the
streamflow gauge station may have contributed to the ranking results. The most downstream
catchment, area 7, ranked first in terms of impact on annual maximum daily flows, and
second in terms of impact on 7Q10 and monthly median daily flows. Catchment area 4
associated with the highest groundwater recharge was ranked first and second for impact on
7Q10 based on a and |3 indices, respectively. Areas characterized by steep topography (area
3) and significant wetlands (area 6) ranked low, some times the lowest, with respect to
impact on the three design flows.
This model study predicted that low and high flows may, respectively, decrease and
increase significantly as a result of urbanization. Land use changes in that part of the
watershed where Bisbing Run and Bulgers Run flow through were predicted to have the
highest impact on the reduction in low flows due to anticipated reduction in the groundwater
recharge rates. The most downstream of the channel network, immediately upstream of the
USGS gauge station, ranked first in terms of impact on monthly median of daily flow and
maximum annual daily flow. Land disturbances in the topographically steep subwatershed, in
which Coolmoor Creek flows, and the area that contains wetlands, through which Cranberry
Creek flows, were predicted to generally have the least impact on the watershed response.
53
-------
8 References
Arnold, J.G., and N. Fohrer. SWAT2000: current capabilities and research opportunities in
applied watershed modeling. Hydrol. Process. 19: 563-572 (2005)
Arnold, J.G., Muttiah, R.S., Srinivasan, R., Allen, P.M. Regional estimation of base flow and
groundwater recharge in the Upper Mississippi river basin. J. Hydrol. 227(1): 21-40
(2000).
Arnold, J.G., and P.M. Allen. Automated methods for estimating base flow and ground water
recharge from streamflow records. J. Am. Water. Resour. As. 35(2): 411-424 (1999).
Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R. Large area hydrologic modeling
and assessment part I: Model development. Journal of the American Water Resources
Association 34(1): 73-89 (1998).
Arnold, J.G. and P.M. Allen. Estimating hydrologic budgets for three Illinois watersheds. J.
Hydrol. 176(1): 57-77 (1996).
Arnold, J.G., Allen, P.M., Muttiah, R., and Bernhardt, G. Automated base flow separation
and recession analysis techniques. Ground Water 33(6): 1010-1018 (1995).
Beven, K. and A. Binley. The future of distributed models: Model calibration and uncertainty
prediction. Hydrol. Process. 6: 279-298 (1992).
Beven, K. and J. Freer. Equifinality, data assimilation, and uncertainty estimation in
mechanistic modeling of complex environmental system using the GLUE methodology. J.
Hydrol. 249: 11-29(2001).
Binley A.M., KJ. Beven, A. Calver, and L.G. Watts. Changing Responses in Hydrology:
Assessing the uncertainty in physically based model predictions. Water Resour. Res. 27(6):
1253-1261 (1991).
Borah, D.K. "Watershed scale non-point source pollution models: mathematical bases." In
Proceedings of the 2002 ASAE Annual International Meeting/ClGR World Congress.,
Chicago, IL, 2002. Paper no. 022091.
Box, G.E.P., and G.M. Jenkins. "Time Series Analysis, Forecasting, and Control" San
Francisco: HoldenDay, 1970.
Carpenter, T.M. and K.P. Georgakakos. Impacts of parametric and radar rainfall uncertainty
on the ensemble streamflow simulations of a distributed hydrologic model. J. Hydrol. 298:
202-221 (2004).
Di Luzio, M., Srinivasan, R., Arnold, J.G., and Neitsch, S.L. ArcView interface for
SWAT2000, User's Guide, TWRI report TR-193. Texas Water Resources Institute, College
Station, Texas, 2002.
Dilks, D.W., R.P. Canale, and P.G. Meier. Development of Bayesian Monte Carlo techniques
for water quality model uncertainty. Ecolog. Model. 62: 149-162 (1992).
Eckhardt, K., and J.G. Arnold. Automated calibration of a distributed catchment model. J.
Hydrol. 251:103-109 (2001).
54
-------
Fohrer, N., Moller, D., and Steiner, N. An interdisciplinary modeling approach to evaluate
the effects of land use change. Phys. Chem. Earth 27:655-662 (2002).
Freer J. and K. Beven. Bayesian estimation of uncertainty in runoff prediction and the value
of data: An application of the GLUE approach. Water Resour. Res. 32(7): 2161-2173
(1996).
Green, W.H., and G.A. Ampt. Studies on soil physics, 1. The flow of air and water through
soils. J. Agr. Sci. 4:11-24 (1911).
Haan, C.T. "Statistical Methods in Hydrology." Second Edition, Iowa State Press, A
Blackwell Publishing Company. 496 pp. 2002.
Hantush, M.M., and L. Kalin. Uncertainty and sensitivity analysis of runoff and sediment
yield in a small agricultural watershed with KINEROS2. Hydro!. Sci. J 50(6): 1151-1171
(2005)
Hossain, F., E.N. Anagnostou, T. Dinku, and M. Borga. Hydrological model sensitivity to
parameter and radar rainfall estimation uncertainty. Hydrol. Process. 18: 3277-3291
(2004).
Hossain, F. and E. N. Anagnostou. Assessment of a probabilistic scheme for flood prediction.
J. Hydrolog. Engrng. 10(2): 141-150 (2005)
Jha, M., Gassman, P.W., Secchi, S., Gu, R., Arnold, J.G., 2003. "Hydrologic simulation of
the Maquoketa River watershed with SWAT." In Proceedings of the AWRA 2003 Spring
Specialty Conference, Kansas City, MO, 2003.
Jayakrishnan, R., Srinivasan, R., Santhi, C., and Arnold, J.G. Advances in the application of
the SWAT model for water resources management. Hydrol. Process. 19:749-762 (2005).
Kalin, L., and M.M. Hantush. Hydrologic Modeling of the Pocono Creek Watershed with
NEXRAD and Rain Gauge Data. J. Hydrolog. Engrng. In press (2006a).
Kalin, L., M.M. Hantush. "Effect of Urbanization on Sustainability of Water Resources in
the Pocono Creek Watershed." In V.P. Singh and Y.J. Xu (ed.) Coastal Hydrology and
Processes, by American Institute of Hydrology, Water Resources Publications, 2006b, 59-
70.
Kalin, L., M.M. Hantush. Evaluation of sediment transport models and comparative
application of two watershed models, EPA/600/R-03/139. Cincinnati, OH: U.S.
Environmental Protection Agency, 2003,
http://www.epa.gov/ORD/NRMRL/Pubs/600r03139/600r03139.pdf
Kim, T-W and J.B. Valdes. Synthetic generation of hydrologic time series based on
nonparametric random generation. J. Hydrol. Engrng. 10(5): 395-404 (2005).
McKay, M.D., Beckman, R.J., and Conover, W. J. A Comparison of Three Methods for
Selecting Values of Input Variables in the Analysis of Output from a Computer Code.
Technometrics 21(2):239-245 (1979).
Lall, U. Recent advances in nonparametric function estimation: Hydrologic applications. Rev.
Geophys. 33(S): 1093-1102 (1995).
55
-------
Muleta, M.K., and J.W. Nicklow. Sensitivity and uncertainty analysis coupled with automatic
calibration for a distributed watershed model. J. Hydrol. 306:127-145 (2005).
Nash, I.E., and Sutcliffe, J.V. River flow forecasting through conceptual models. Part 1: A
discussion of principles. J. Hydrol. 10:282-290 (1970).
Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., and King, K.W. Soil and Water
assessment Tool Theoretical Documentation, version 2000. TWRI report TR-191, Texas
Water Resources Institute, College Station, Texas, 2002a.
Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., and King, K.W. Soil and Water
Assessment Tool User's Manual, version 2000. TWRI report TR-192, Texas Water
Resources Institute, College Station, Texas, 2002b.
NRC. "Assessing the TMDL approach to water quality management'' Committee to Access
the Scientific Basis of the Total Maximum Daily Load Approach to Water Pollution
Reduction, Water Science and Technology Board, Division on Earth and Life Studies,
National Research Council, Washington, D.C, 2001.
Pebesma, E.J., P. Switzer, and K. League. Error analysis for the evaluation of model
performance: rainfall-runoff event time series data. Hydrol. Process. 19: 1529-1548 (2005).
Pocono Creek Pilot Study Draft Technical Report. Geology, geomorphology, geohydrology,
and surface water hydrology of the Pocono Creek Basin,
(Jan. 30, 2006). Delaware River Basin
Commission (DRBC), 2001.
Rallison, R.E., and Miller, N. "Past, present and future SCS runoff procedure." In Rainfall
runoff relationship. V.P. Singh, eds., Water Resources Publication, Littleton, Colo., 1981,
353-364.
Saghafian, B. and Khosroshahi, M. Unit response approach for priority determination of
flood source areas. J. Hydro. Engrng. 10(4): 270-277 (2005).
Sal eh, A. and B. Du. "Application of SWAT and HSPF within BASINS program for the
Upper North Bosque River watershed." In Proceedings of the 2002 ASAE Annual
International Meeting, 2002.
Saltelli, A., K. Chan, and E.M. Scott. 2000. "Sensitivity Analysis" John Wiley & Sons, Ltd,
2000, 475 pp.
Santhi, C., Arnold, J.G, Williams, J.R., Dugas, W.A., Srinivasan, R., and Hauck, L.M.
Validation of the SWAT model on a large river basin with point and nonpoint sources. J.
Am. Water. Resour. As. 37(5): 1169-1188 (2001).
Sharply, A.N., and J.R. Williams. Eds. EPIC-Erosion Productivity Impact Calculator, 1.
Model Documentation, U.S. Department of Agriculture, Agricultural Research Service,
Tech. Bull. 1768, 1990.
Shumway, R.H. "Applied Statistical Time Series Analysis" Prentice Hall, 1988, 379 pp.
Silverman, B.W. "Density Estimation for Statistics and Data Analysis" Chapman and Hall,
New York, 1986.
Smakhtin, V.U. Low Flow Hydrology: A Review. J. Hydrol. 240:147-186 (2001).
56
-------
Sophocleous, M., and S.M. Perkins. Methodology and application of combined watershed
and ground-water models in Kansas. JHydrol. 236:185-201 (2000).
Sorooshian, S. and J.A. Dracup. Stochastic parameter estimation procedures for hydrologic
rainfall-runoff models: Correlated and Heteroscedastic error cases. Water Resour. Res.
16(2): 430-442 (1980).
Sorooshian, S., Gupta, V.K., and Fulton, J.L. Evaluation of maximum likelihood parameter
estimation techniques for conceptual rainfall-runoff models: Influence of calibration data
variability and length on model credibility. Water Resour. Res. 19(1):251-259 (1983).
Spear, R.C. and G.M. Hornberger. Eutrophication in peel inlet-II. Identification of critical
uncertainties via generalized sensitivity analysis. Water Res. 14: 43-49 (1980).
Spear, R.C., T.M. Grieb, and N. Shang. Parameter uncertainty and interaction in complex
environmental models. Water Resour. Res. 30(11): 3159-3169 (1994).
Srinivasan, R., Ramanarayanan, T.S., Arnold, J.G., Bednarz, S.T., Large area hydrologic
modeling and assessment part II: model application. Journal of the American Water
Resources Association 34(1): 91-101 (1998).
Tripathi, M.P., Panda, R.K., Raghuwanshi, N.S., and Singh, R. Hydrological modeling of a
small watershed using generated rainfall in the soil and water assessment tool model.
Hydrol. Process. 18:1811-1821 (2004).
US Department of Commerce, Census Bureau, Statistical Abstract of the United States, 2001.
URL: http://www.census.gov/prod/2005pubs/06statab/pop.pdf, 2005.
USDA Soil Conservation Service. "National Engineering Handbook Section 4 Hydrology."
Chapters 4-10, 1972.
Van der Perk, M. and M.F.P. Bierkens. The identifiability of parameters in a water quality
model of the Biebrza, Poland. J. Hydrol. 200: 307-322 (1997).
Van Liew, M.W., Arnold, J.G., Gardbercht, J.D. Hydrologic simulation on agricultural
watersheds: choosing between two models. Trans. ASAE 46(6): 1539-1551 (2003).
Williams, J.R. Flood routing with variable travel time or variable storage coefficients. Trans.
ASAE 12(1): 100-103 (1969).
57
------- |