United States
Environmental Protection
Agency
•^^
for Selection and Spatial Allocation of
Non-Point Source Pollution Control Practices
fi
Mazdak Arabi, Rao S, Govhdaraju, Mohamed M, Hantush
-------
EPA/600/R-08/036
January 2007
WATERSHED MANAGEMENT TOOL FOR SELECTION
AND SPATIAL ALLOCATION OF NONPOINT SOURCE
POLLUTION CONTROL PRACTICES
Prepared By
Mazdak Arabi
Department of Agricultural and Biological Engineering
Purdue University
West Lafayette, Indiana 47907
Rao S. Govindaraju
School of Civil Engineering
Purdue University
West Lafayette, Indiana 47907
Contract Number: 4C-R330-NAEX
Project Officer
Mohamed M. Hantush
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. This report has been subjected to the Agency's peer and
administrative review and has been approved for publication as an EPA document. 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 water quality
in eutrophic water bodies and did not involve collection and analysis of environmental
measurements.
11
-------
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
in
-------
ABSTRACT
Distributed-parameter watershed models are often utilized for evaluating the effectiveness of
sediment and nutrient abatement strategies through the traditional (calibrate^- validate^-
predict} approach. The applicability of the method is limited due to modeling
approximations. In this study, a computational method is presented to determine the
significance of modeling uncertainties in assessing the effectiveness of best management
practice (BMPs) in two small watersheds in Northeastern Indiana with the Soil and Water
Assessment Tool (SWAT). The uncertainty analysis aims at (i) identifying the hydrologic
and water quality processes that control the fate and transport of sediments and nutrients
within watersheds, and (ii) establishing uncertainty bounds for model simulations as well as
estimated effectiveness of BMPs. The SWAT model is integrated with a Monte-Carlo based
methodology for addressing model uncertainties. The results suggested that fluvial processes
within the channel network of the study watersheds control sediment yields at the outlets, and
thus, BMPs that influence channel degradation or deposition are the more effective sediment
control strategies. Conversely, implementation of BMPs that reduce nitrogen loadings from
uplands areas such as parallel terraces and field borders appeared to be more crucial in
reducing total N yield at the outlets. The uncertainty analysis also revealed that the BMPs
implemented in the Dreisbach watershed reduced sediment, total P, and total N yields by
nearly 57%, 33%, and 31%, respectively. Finally, a genetic algorithm (GA)-based
optimization methodology is developed for selection and placement of BMPs within
watersheds. The economic return of the selected BMPs through the optimization model was
nearly three fold in comparison to random selection and placement of the BMPs.
IV
-------
TABLE OF CONTENTS
Notice ii
Foreword iii
Abstract iv
List of Figures viii
List of Tables ix
Acknowledgements xi
SECTION 1. Introduction 1
1.1. Rationale and Background 1
1.2. Objectives 5
1.3. Impacts of the Study 5
1.4. Structure of the Report 6
SECTION 2. Case Study Watersheds and the Watershed Model 7
2.1. Case Study Watersheds 7
2.2. Choice of the Watershed Model: SWAT 8
2.2.1. Model calibration 11
2.2.2. Representation of BMPs 12
SECTION 3. Probabilistic Approach for Analysis of Uncertainty in the Evaluation of
Watershed Management Practices 14
3.1. Abstract 14
3.2. Introduction 15
3.3. Theoretical Considerations 17
3.3.1. Morris's One-at A Time sensitivity analysis 17
3.3.2. Generalized Likelihood Uncertainty Estimation (GLUE) 17
3.4. Methodology 18
v
-------
3.4.1. OAT sensitivity analysis 19
3.4.2. Range adjustment 20
3.4.3. Uncertainty computations 21
3.5. Results 22
3.5.1. Range adjustment results 22
3.5.2. Uncertainty analysis results 26
3.6. Discussion 31
3.7. Conclusions 32
SECTION 4. Sensitivity Analysis of Sediment and Nutrient Processes with a Watershed
Model 33
4.1. Abstract 33
4.2. Introduction 33
4.3. Theoretical Considerations 36
4.3.1. Regionalized Sensitivity Analysis (RSA) 36
4.3.2. Tree Structured Density Estimation (TSDE) 36
4.4. Methodology 38
4.4.1. Behavior Definition 39
4.5. Results 40
4.6. Discussion 43
4.7. Conclusions 44
SECTION 5. Cost-Effective Allocation of Watershed Management Practices Using a Genetic
Algorithm 54
5.1. Abstract 54
5.2. Introduction 54
5.3. Theoretical Considerations 58
5.3.1. Genetic Algorithm (GA) 58
5.4. Methodology 59
5.4.1. NFS Model 60
5.4.2. BMP Representation 61
5.4.3. Economic Component 61
5.4.4. Optimization Component 64
5.4.5. Optimization Cases 70
VI
-------
5.5. Results 71
5.6. Discussion 77
5.7. Conclusions 78
SECTION 6. Conclusions 80
Bibliography 83
vn
-------
LIST OF FIGURES
Figure Page
2.1 Case Study Watersheds 8
3.1 Spider plots for the most sensitive input factors 24
3.2 Cumulative GLUE-likelihood for output variables simulated at the outlet of the
Dreisbach watershed 27
3.3 Cumulative GLUE-likelihood for output variables simulated at the outlet of the Smith
Fry watershed 28
4.1 Computational framework for sensitivity analysis of sediment and nutrient processes....39
4.2 Posterior distributions of behavior and non-behavior sets of input factors along with
their prior distribution (uniform distribution) for sediment at the outlet of Dreisbach
and Smith Fry watersheds 46
4.3 Posterior distributions of behavior and non-behavior sets of input factors along with
their prior distribution (uniform distribution) for total N at the outlet of Dreisbach and
Smith Fry watersheds 47
4.4 TSDE diagram for sediment at the outlet of Dreisbach watershed 48
4.5 TSDE diagram for sediment at the outlet of Smith Fry watershed 48
4.6 TSDE diagram for total N at the outlet of Dreisbach watershed 49
4.7 TSDE diagram for total Nat the outlet of Smith Fry watershed 49
5.1 Schematic of the optimization procedure 60
5.2 Schematic of an optimization string 65
5.3 Analysis of sensitivity of GA operating parameters 71
5.4 (a) Case3 optimization outputs for the Dreisbach watershed over 2000-2009 period.
(b) Case 4 optimization outputs for the Dreisbach watershed over 2000-2009 period.
(c) Case 4 optimization outputs for the Smith Fry watershedover2000-2009period 74
5.5 Spatial allocation of BMPs from: (a) case 3 computations for Dreisbach watershed;
(b) case 4 computations for Dreisbach; (c) case 4 computations for smith Fry
watershed 76
5.6 Tradeoff between economic criterion (cost) and environmental criteria 77
Vlll
-------
LIST OF TABLES
Table Page
2.1 BMPs installed in the Dreisbach and Smith Fry watersheds 8
2.2 Listing of SWAT parameters 10
2.3 Representation of BMPs with SWAT 13
3.1 Top five sensitive SWAT parameters for streamflow, sediment, total P, and total N
computations 23
3.2 Adjusted parameter ranges for the study watersheds 25
3.3 Summary of the results from analysis of uncertainty 29
4.1 Regionalized Sensitivity Analysis ranking of input factors for sediment based on
observed-past behavior definition at the outlet of Dreisbach watershed 50
4.2 Regionalized Sensitivity Analysis ranking of input factors for sediment based on
observed-past behavior definition at the outlet of Smith Fry watershed 50
4.3 Regionalized Sensitivity Analysis ranking of input factors for total N based on
observed-past behavior definition at the outlet of Dreisbach watershed 51
4.4 Regionalized Sensitivity Analysis ranking of input factors for total N based on
observed-past behavior definition at the outlet of Smith Fry watershed 51
4.5 Summary of TSDE terminal nodes for sediment related parameters based on
observed-past behavior definition at the outlet of Dreisbach watershed 52
4.6 Summary of TSDE terminal nodes for sediment related parameters based on
observed-past behavior definition at the outlet of Smith Fry watershed 52
4.7 Summary of TSDE terminal nodes for total N related parameters based on observed-
past behavior definition at the outlet of Dreisbach watershed 53
4.8 Summary of TSDE terminal nodes for total N related parameters based on observed-
past behavior definition at the outlet of Smith Fry watershed 53
5.1 Unit cost of establishment (c), maintenance rate (rm), and design life (td) of BMPs in
the study along with number (nbmp) and quantity (abmp) of each type in the study
watersheds 62
IX
-------
5.2 Monetary benefits of unit reduction of sediment and nutrient loads in 2000 dollars 64
5.3 Various combinations of GA operating parameters for sensitivity analysis 69
5.4 Description of cases compared in this study 70
5.5 Results for base-case and targeting cases. Mean and maximum value of monthly
sediment, phosphorus, and nitrogen loads were estimated from SWAT simulations ....72
5.6 Results for optimization cases 73
x
-------
ACKNOWLEDGEMENTS
The U.S. Environmental Protection Agency through its Office of Research and Development
funded the research described here. This support is acknowledged. The report benefited from
the constructive comments of Dr. M. M. Hantush. The report was reviewed by Drs.
Zhonglong Zhang and M. Dougherty, and their comments were utilized in preparing this
final version. Their help is acknowledged.
XI
-------
SECTION 1. INTRODUCTION
1.1. Rationale and Background
Soil and water conservation practices are widely accepted as effective control measures for
agricultural nonpoint sources of sediments and nutrients. The 2002 Farm Bill provided up to
$13 billion for conservation programs aimed at protecting water quality from agricultural
nonpoint source (NFS) pollution (USDA, 2003). In addition, under the Clean Water Act
Section 319 Nonpoint Source National Monitoring Program and wetland protection
programs, the EPA supports programs to reduce the negative impacts of runoff from
agricultural, urban, and industrialized areas. Similarly, the Natural Resources Conservation
Service (NRCS) provides hundreds of millions of dollars in federal funds to support
agricultural conservation practices in an effort to reduce the movement of pollutants into our
waterways. Success of such programs, however, is contingent upon availability of efficient
watershed-scale planning tools.
Watershed management plans are composed of individual structural and/or cultural soil and
water conservation units that are often referred to as best management practices (BMPs).
BMPs can substantially reduce transport of contaminants to surface water (Mostaghimi et al.,
1997; Arabi et al., 2004), but, their implementation bears additional costs to stakeholders and
watershed managers. Assessment of watershed management plans is challenged by
complexities in incorporation of conflicting environmental, economic, and institutional
criteria. Environmental assessments in watersheds hinge on resolving social benefits such as
achieving the goal of swimable and fishable water bodies under the US EPA's total
maximum daily load (TMDL) agenda. While BMPs facilitate achievement of such targets,
their establishment bears additional cost for watershed management and/or agricultural
producers. Since management practices are usually implemented under a limited budget,
costs associated with unnecessary/redundant management actions may jeopardize
attainability of designated water quality goals. Identifying optimal combinations of
-------
watershed management practices requires systematic approaches that allow decision makers
to quickly assess trade-offs among environmental and economic criteria.
Implementation of BMPs is rarely followed by long-term monitoring and data collection
efforts to assess their effectiveness in meeting their intended goals. Long-term data on flow
and water quality within watersheds, before and after placement of BMPs, is not generally
available. Therefore, cost-effective evaluation of BMPs (especially new ones that have had
little or no history of use) needs to be conducted with the help of simulation models. The
assessment of these economic and environmental metrics requires an integrated approach
often involving social objectives and competing interests from various stakeholder groups.
Arabi et al. (2004) developed a processed based methodology, using the Soil and Water
Assessment Tool (SWAT), to numerically represent the impacts of four structural BMPs on
hydrologic/water quality processes in two small agricultural watersheds in Indiana. Various
processes that were considered included: infiltration; surface runoff (peak and volume);
upland erosion (sheet and rill erosion); gully and channel erosion; nutrient loadings from
upland areas; and within-channel processes. Based on the function of a BMP, specific SWAT
parameters that represent the impacted hydrologic/water quality processes were altered.
While similar approaches can be followed for other BMPs, an important question for
evaluation of BMPs at the watershed scale arises: how can the tools of modeling and
computational analysis be applied to identify the best configuration (type and location) of
BMPs at the watershed scale!
As stated earlier, this question calls for an integrated assessment involving process-based
models, economics, and social considerations. Consequently, several models have to be
utilized in conjunction to represent the various processes that will have a bearing on the
modeling targets and the decision making process. Computational techniques can assist with
quantification of modeling uncertainties, critical natural processes and key management
actions, and optimal location of BMPs.
Uncertainty Analysis
Any modeling process will necessarily entail reducible and inherent uncertainty from data,
model abstractions, and natural heterogeneity of watersheds. The National Research Council
-------
report "Assessing the TMDL approach to water quality management; NRC, 2001) concluded
that modeling uncertainty should be rigorously and explicitly addressed in development and
application of models for environmental management, especially when stakeholders are
affected by the decisions contingent upon model-supported analyses. The EPA's guideline
for TMDL development recommends accounting for the uncertainties embedded in model
estimates by applying a margin of safety (MOS), i.e., TMDL=^WLA+^LA+MOS where
7MDZ=maximum pollutant load a water body can receive and still maintain water quality
standards, ^JVLA=pomt source waste load allocation, and ^LA= nonpoint source load
allocation. Lack of systematic and efficient approaches for evaluation of uncertainties
associated with cost-effectiveness of watershed plans has led to abandonment of explicit
computation of the MOS (Dilks and Freedman, 2004).
The uncertainty estimates (i.e., MOS values) associated with absolute estimates of design
quantities of agricultural NFS pollution tend to be very high because of data sparsity and
model limitations (Osidele et al., 2003; Benaman and Shoemaker, 2004). Although the
literature is replete with sensitivity analysis and uncertainty analysis methods, implications of
uncertainty associated with model predictions have not been widely endorsed in the decision
making process mainly as a result of large uncertainty estimates. The magnitude of
uncertainty itself is a key factor in its acceptance as the cost of implementation of
management actions such as the TMDL program may significantly increase with larger
uncertainty estimates (Dilks and Freedman, 2004). The basic thesis of this report is that if the
goal of an integrated modeling study is to examine the impact of management scenarios on
water quality of a study area, it may be neither practical nor necessary to incorporate large
uncertainty of absolute predictions in the decision making process. It would be more
feasible, and desirable, to communicate uncertainty of estimated effectiveness of management
scenarios rather than uncertainty of absolute predictions.
Identification of critical natural processes for selection of management actions
The current US EPA's "draft handbook for developing watershed plans to restore and protect
our waters" (US EPA, 2005a) recommends implementation of management practices in
portions of the watershed that are believed to contribute intensively to NPS pollution.
Locating critical areas, however, is complicated because contaminants are carried along with
-------
the flow, and water movement over a watershed tends to be fairly dynamic, behaving in a
nonlinear fashion. Potential interactions among hydrologic, fluvial, and nutrient processes
and NFS control measures should be identified and included in designing watershed
management scenarios. Nearly a hundred BMPs for NFS pollution control are recommended
by the USD A Natural Resources Conservation Service (USD A NRCS, 2006). Evaluation of
impacts of all these practices, even within small watersheds, is infeasible. However, it is
likely that only a few management actions largely control fate and transport of contaminants
within a given watershed system (Arabi et al., 2006a). We hypothesize that the analysis of
uncertainty of model simulations could highlight critical processes and key management
actions that are key to control of NFS pollution in a given watershed system.
Identification of optimal spatial configuration of management actions
Management of limited and fragile environmental resources in the face of conflicting
interests is a challenging problem. An appropriate balance needs to be struck between
protecting the environment and allowing users to access the natural resource base and the
environment to assure their livelihoods. For example, farmers need to be able to pursue
avenues that increase agricultural productivity without creating adverse environmental
consequences. This includes safeguarding their plants from pests, augmenting soil fertility to
raise productivity to levels required in today's modern agriculture, and accessing water.
Systematic methods are required to assist policy makers in selecting strategies that lead to
environmentally friendly solutions along with being cost-effective. To make informed
decisions, policy makers should have the tools necessary for evaluation of alternatives and to
assess trade-offs between environmental and socioeconomic criteria.
Research to date indicates the promise of heuristic optimization for cost-effective allocation
of watershed management practices (Srivastava et al., 2002; Veith et al., 2004; Muleta and
Nicklow, 2005). Unlike gradient-based approaches, heuristic techniques do not require
linearity, continuity, or differentiability either for objective/constraint functions or for input
parameters. Thus, they are well-suited for cost-effective allocation of watershed management
plans. We hypothesize that heuristic optimization techniques could facilitate spatial
allocation of management actions in watersheds in a more cost-effective fashion than
random and even more commonly used cost-sharing and targeting strategies.
-------
1.2. Objectives
The overall goal of this hypothesis-driven study is to develop a computational framework for
selection and placement of BMPs in a cost-effective manner. Ultimately, the hope is that this
tool will help watershed management programs, such as the TMDL program, attain the
desired water quality goals for impaired water bodies. This goal is accomplished through the
following specific objectives:
1. Development of an uncertainty analysis method for establishing uncertainty bounds for
the estimated effectiveness of BMPs for sediment and nutrient control.
2. Development of a computational framework for identification of natural fluvial processes
and BMPs that control fate and transport of sediments and nutrients at the watershed
scale.
3. Development of an optimization tool for cost-effective allocation of BMPs within
watersheds.
1.3. Impacts of the Study
The impacts of the present work in management of watershed systems are three fold. First,
the computational analysis presented here will increase the capabilities of watershed
managers to quantify and integrate modeling uncertainties in the decision making process.
The tool is likely to promote stakeholder involvement through cost- a language that they
understand and are able to identify with. The additional cost of implementation of watershed
management plans due to incorporation of modeling uncertainties is key to adoption of
scientific uncertainty estimates in the planning process. The present work presents an
effective framework for analysis of uncertainties and its communication to stakeholders.
Second, the analysis will allow watershed managers to asses what it takes to reduce pollutant
loads/concentrations to meet water quality standards. For example, the computational
analysis can be used to identify key management actions to reduce average annual nitrate
concentration below 10 ppm, the EPA's drinking water standard. Third, the optimization
framework presented herein could greatly increase the water quality benefits of dollars spent
on conservation plans. Alternatively, the analysis could be used to attain the same level of
-------
water quality benefits currently received from watershed management plans for substantially
less costs.
1.4. Structure of the Report
The remainder of this report is organized into five sections. Section 2 provides a description
of the study watersheds and the choice of the watershed model. The hypotheses and
objectives of this study will be tested in two small agricultural watersheds in northeastern
Indiana within the Maumee River basin located in the Midwestern portion of the U.S. A
review of the choice of the watershed model, Soil and Water Assessment Tool (SWAT) will
also be presented. In Section 3, a 3-step computational analysis based on Latin Hypercube
Sampling, the Morris's One-at-A-Time sensitivity analysis and the Generalized Likelihood
Uncertainty Estimation (GLUE) is developed and demonstrated for establishing uncertainty
bounds associated with evaluation of BMPs. Application of multivariate multiobjective
sensitivity analysis including Regionalized Sensitivity Analysis (RSA) and Tree-Structured
Density Estimation (TSDE) for identification of critical natural processes and key
management actions for sediment and nutrient control is presented in Section 4. This is
followed in Section 5 by development and demonstration of a genetic algorithm (GA) search
engine for identification of best location of management actions/BMPs at the watershed
scale. Finally, major conclusions of the study and recommendations for future research are
discussed in Section 6.
The report is primarily based on three journal papers: Sections is derived from Arabi et al.
(2007a), Section 4 is derived from Arabi et al. (2007b), and Section 5 is derived from Arabi
et al. (2006a).
-------
SECTION 2. CASE STUDY WATERSHEDS AND THE WATERSHED MODEL
2.1. Case Study Watersheds
The utility of the optimization framework was examined for two subwatersheds in the Black
Creek basin. The Black Creek watershed (Figure 2.1) located in northeast Indiana is a typical
watershed in the upper Maumee River basin in the Midwestern portion of the United States.
In mid 1970's and early 1980's, several BMPs were implemented in the watershed and
detailed water quality monitoring was carried out at various locations within the watershed to
examine short-term water quality impacts of soil and water conservation techniques (Lake
and Morrison, 1977; Lake and Morison, 1978). Data collected from automated samplers at
the outlets of Dreisbach (6.23 km2) and Smith Fry (7.3 km2) within the Black Creek
watershed were the most complete and were used in this study. Black Creek watershed has
been listed as an impaired water body in Indiana with nutrients, algal growth, and impaired
biotic communities as major concerns (EPA, 2005 b).
The dominant hydrological soil group in the study watersheds is type C. Land use in
Dreisbach is mostly pasture in the upper portion, while row crops are wide spread in the
remainder of the watershed. Smith Fry is a predominantly agricultural watershed. Field
borders, parallel terraces, grade stabilization structures, and grassed waterways were the
structural BMPs installed in the watershed by targeting (see Figure 2.1). Table 2.1
summarizes the number and area under influence of each type of BMP installed in the study
watersheds as shown in Figure 2.1. Available data, land use distribution, and other
information for the watersheds can be obtained from Lake and Morrison (1977), Lake and
Morison (1978), and Morrison and Lake (1983), and have been more recently described in
Arabi et al. (2004).
-------
/Lake Erie
. y--
S£T;
* •* Maumee River
X N
\
\ Black Creek
watershed
State of
Indiana
0 100 200
400
Kilometers
# Field Border
$ Parallel Terrace
N Grade Stabilization Structure
N Grassed Waterway
Kilometers
Figure 2.1 Case Study Watersheds.
Table 2.1 BMPs installed in the Dreisbach and Smith Fry watersheds
BMP
Field Border
Parallel Terrace
Grassed Waterway
Grade Stabilization Structure
Dreisbach
,T , Length/Area
Number ? . ,
(unit)
7
4
5
10
2600 (m)
2130 (m)
3. 50 (ha)
Smith Fry
,T , Length/Area
Number ? . ,
(unit)
1
2
1
2
1800(m)
480 (m)
0.95 (ha)
2.2. Choice of the Watershed Model: SWAT
Soil and Water Assessment Tool (SWAT) (Arnold et al., 1998; Arnold and Fohrer, 2005) is a
process-based distributed-parameter simulation model, operating on a daily time step. The
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. SWAT uses readily available inputs and has the capability of routing
runoff and chemicals through streams and reservoirs, and allows for addition of flows and
inclusion of measured data from point sources. Moreover, SWAT has the capability to
evaluate the relative effects of different management scenarios on water quality, sediment,
and agricultural chemical yield in large, ungaged basins. Major components of the model
include weather, surface runoff, return flow, percolation, evapotranspiration (ET),
transmission losses, pond and reservoir storage, crop growth and irrigation, groundwater
flow, reach routing, nutrient and pesticide loads, and water transfer. Table 2.2 provides a
listing of important SWAT input parameters corresponding to the above-mentioned
components. The lower and upper parameter bounds (LB and UP) were obtained from
SWAT theoretical documentation (Neitsch et al., 2002) and similar studies in the literature
(Lenhart et al., 2002; Benaman and Shoemaker, 2004; van Griensven et al., 2006).
For simulation purposes, SWAT partitions the watershed into subunits including subbasins,
reach/main channel segments, impoundments on main channel network, and point sources.
Subbasins are divided into hydrologic response units (HRUs) that are portions of subbasins
with unique land use, management, and soil attributes. SWAT uses a modification of the SCS
curve number method (USDA Soil Conservation Service, 1972) or Green and Ampt
infiltration method (Green and Ampt, 1911) to compute surface runoff volume for each
HRU. Peak runoff rate is estimated using a modification of the Rational Method. Daily or
sub-daily rainfall data is used for calculations. Flow is routed through the channel using a
variable storage coefficient method developed by Williams (1969) or the Muskingum routing
method.
Erosion and sediment yield are estimated for each FtRU with the Modified Universal Soil
Loss Equation (MUSLE) (Williams, 1975). Sediment deposition and degradation are the two
dominant channel processes that affect sediment yield at the outlet of the watershed. Whether
channel deposition or channel degradation occurs depends on the sediment loads from upland
areas and transport capacity of the channel network. If sediment load in a channel segment is
larger than its sediment transport capacity, channel deposition will be the dominant process.
Otherwise, channel degradation (i.e. channel erosion) occurs over the channel segment.
-------
Table 2.2 Listing of SWAT parameters: LB and UB refer to lower and upper bounds of
parameter vector, parameters identified by were altered as a percentage of the default
value, parameters identified by " were considered in sensitivity analysis but not in
uncertainty analysis.
SWAT Symbol Description
CN2'
SOL_AWC
ESCO
OV N
SLOPE"' "
SLSUBBSN
GWQMN
REVAPMN
GW_REVAP
GW_DELAY
ALPHA_BF
SURLAG
SFTMP
CH_Sl""
CH_N1
CH_K1
CH_S2"'"
CH_N2
CH_K2
CH_EROD
CH_COV
PRF
SPCON
SPEXP
BIOMIX
BIOMIN
USLE_P
USLE_C"
USLE_K
ORGP_AG
ORGP_PAST
LABP_AG
LABP_PAST
ORGN_AG
ORGN_PAST
SOLN_AG
SOLN_PAST
SCS runoff curve number
available soil water capacity
soil evaporation compensation factor
Manning's "n" value for overland flow
average slope steepness
average slope length
minimum threshold depth of water in the shallow aquifer for return
minimum threshold depth of water in the shallow aquifer for "revap"
groundwater "revap" coefficient
groundwater delay
baseflow alpha factor for recession constant
surface runoff lag time
snowfall temperature
average slope for tributary channels
Manning's "«" value for tributary channels
effective hydraulic conductivity in tributary channels
average slope for the main channels
Manning's "«" value for the main channel
effective hydraulic conductivity in the main channel
channel credibility factor
channel cover factor
peak rate adjustment factor for in-stream channel routing
Linear coefficient for in-stream channel routing
exponent coefficient for in-stream channel routing
biological mixing efficiency
minimum plant biomass for grazing
USLE equation support practice factor
minimum value of USLE equation cover factor
USLE equation soil credibility factor
initial organic P in soils for agriculture land use
initial organic P in soils for pasture land use
initial soluble P in soils for agriculture land use
initial soluble P in soils for pasture land use
initial organic N in soils for agriculture land use
initial organic N in soils for pasture land use
initial NO3 in soils for agriculture land use
initial NO3 in soils for pasture land use
Units
%
mm/mm
m/m
m
mm
mm
days
days
oC
m/m
Range
LB UB
-15
0.01
0
0.1
0
0.15
0
0
0.02
0
0
1
-5
0
15
0.4
1
0.3
0.6
1.2
5000
500
2
500
1
12
5
10
0.008 0.065
mm/hr
m/m
mm/hr
cm/hr/Pa
kg/ha
%
mg/kg
mg/kg
mg/kg
mg/kg
mg/kg
mg/kg
mg/kg
mg/kg
0
0
0.01
0
0
0
0
0
1
0.01
0
0.1
0.001
0.01
1
1
1
1
1
1
0.1
0.1
150
10
0.3
150
0.6
0.6
2
0.001
1.5
1
1
1
0.5
0.65
1000
500
100
50
10000
5000
5
3
10
-------
Movement and transformation of several forms of nitrogen and phosphorus over the
watershed are accounted for within the SWAT model. Nutrients are introduced into the main
channel through surface runoff and lateral subsurface flow, and transported downstream with
channel flow. Major phosphorous sources in mineral soil include organic phosphorus
available in humus and mineral phosphorus that is not soluble. Phosphorus may be added to
the soil in the form of fertilizer, manure, and residue application. Surface runoff is deemed as
the major mechanism of phosphorus removal from a field. Unlike phosphorus that has low
solubility, nitrogen is highly mobile. Major nitrogen sources in mineral soil include organic
nitrogen available in humus, mineral nitrogen in soil colloids, and mineral nitrogen in
solution. Nitrogen may be added to the soil in the form of fertilizer, manure, or residue
application. Plant uptake, denitrification, volatilization, leaching, and soil erosion are the
major mechanisms of nitrogen removal from a field.
It is worthwhile mentioning that the computational analyses proposed in this study do not
necessarily require using the SWAT model as the NFS component. We have chosen SWAT
because the case studies that will be performed will deal with sediments and nutrients.
Soundness of mathematical representation of sediment, nutrient, and pesticide processes in
SWAT has been validated in previous research (Santhi et al., 2001; Kirsch et al., 2002; Arabi
et al, 2004; Vazquez-Amabile and Engel, 2006). Moreover, SWAT has been linked with
optimization routines for automated calibration of the model (van Griensven and Bauwens,
2003) and for evaluation of efficiency of nonpoint source pollution regulatory programs
(Whittaker et al., 2003). The optimization framework herein can be easily integrated with
other hydrologic/water quality models to develop management plans for other types of
pollutants.
2.2.1. Model calibration
Arabi et al. (2004) utilized a procedure adopted from Santhi et al. (2001) to calibrate and
validate the SWAT model for the Dreisbach and Smith Fry watersheds. The results from this
previous study indicated acceptable agreement between observed and simulated baseflow,
streamflow, sediment yield, and total P and total N loads at the outlet of the study
11
-------
watersheds. In this study, the adjusted value of SWAT parameters from the manual
calibration in Arabi et al. (2004) were used as default values.
2.2.2. Representation of BMPs
For this study, a method presented by Arabi et al. (2004) was utilized to evaluate the water
quality impacts of grassed waterways, grade stabilization structures, field borders and
parallel terraces. The method was developed based on published literature pertaining to BMP
simulation in hydrological models and considering the hydrologic and water quality
processes simulated in SWAT. Based on the function of the BMPs and hydrologic and water
quality processes that are directly modified by their implementation, corresponding SWAT
parameters were selected and altered. The impacts of the BMPs on other hydrologic/water
quality processes were indirectly accounted for through their functional representation within
the SWAT model. For example, implementation of grassed waterways prevents gully erosion
due to concentrated flow. Thus, channel cover and erodibility factors were set to zero for
channel segments with grassed waterways. Also, grassed waterways increase deposition of
sorbed pollutant loads in the channel network ((Fiener and Auerswald, 2003). This impact is
captured by increasing the channel Manning's number to 0.3 (Fiener and Auerswald, 2006).
For parallel terraces, Wischmeier and Smith (1978) recommended a USLE practice factor of
0.1 for terraces on slopes less than 15% with graded channels sod outlets. Table 2.3
summarizes SWAT parameters and their corresponding values for representation of the
BMPs. More information on BMP representation procedure with SWAT model can be
obtained from Arabi et al. (2004), and Arabi et al. (2006b).
It is worthwhile to mention that the numerical representation of BMP performances as
presented in Table 2.3 bears some degree of uncertainty. For example, Fiener and Auerswald
(2006) assumed CH_N2 ranges between 0.3-0.4 s m"1/3 over the year for grassed waterways
in case of dense grasses and herbs under non-submerged conditions. Assuming a fixed value
for the representative parameters adds additional uncertainty in the BMP evaluation
methodology. This uncertainty will perhaps be alleviated once precise numerical
representation procedures are developed and validated with experimental data. We have
12
-------
adopted the methodology provided by Arabi et al. (2004) without explicitly incorporating the
uncertainty of BMP representation procedure.
Table 2.3 Representation of BMPs with SWAT
BMP
Function
Representing SWAT Parameter
Parameter
(SWAT input file)
Range
value when BMP
implemented
Field Border increase sediment trapping
FILTERW
(.hru)
0-5 (m)
5(m)
Parallel
Terrace
reduce overland flow
reduce sheet erosion
reduce slope length
CN2
(•mgt)
USLEP
(•mgt)
SLSUBBSN
(.hru)
0-100
0-1
10-150
0.1
(terraced)
Grassed
Waterway
increase channel cover
reduce channel erodibility
increase channel roughness
CH_COV
(.rch)
CH_EROD
(.rch)
CH_N2
(.rch)
0-1
0-1
0-0.3
0.0
(fully protected)
0.0
(non-erosive)
0.24
Grade reduce gully erosion
Stabilization
Structure reduce slope steepness
CH_EROD
(.rch)
CHJ52
(.rch)
0-1
0.0
(non-erosive)
- Estimated based on land use and hydrologic soil group of the HRU where it is installed for
terraced condition.
- Estimated for each parallel terrace based on its features and SWAT assigned overland
slope of the HRU where it is installed:
SLSUBBSN = (A xS+B) x JOO/S
where S is average slope of the HRU; ,4=0.21 and B=0.9 (ASAE, 2003).
- Estimated for each grade stabilization structure based on its features and SWAT assigned
slope, CH_S20id, and length of the channel segment where it is installed:
CH_S2 new = CH S2 oU -D/CH L2
where D is height of the structure (1.2 m) and CH L2 is length of the channel segment.
13
-------
SECTION 3. PROBABILISTIC APPROACH FOR ANALYSIS OF UNCERTAINTY IN THE
EVALUATION OF WATERSHED MANAGEMENT PRACTICES
3.1. Abstract
A computational framework is presented for analyzing the uncertainty in model estimates of
water quality benefits of best management practices (BMPs) in two small (<10 km2)
watersheds in Indiana. The analysis specifically recognizes the significance of the difference
between the magnitude of uncertainty associated with absolute hydrologic and water quality
predictions, and uncertainty in estimated benefits of BMPs. The Soil and Water Assessment
Tool (SWAT) is integrated with Monte Carlo-based simulations, aiming at (1) adjusting the
suggested range of model parameters to more realistic site-specific ranges based on observed
data, and (2) computing a scaled distribution function to assess the effectiveness of BMPs. A
three-step procedure based on Latin Hypercube Sampling (LHS), the Morris's One-factor-
At-a-Time (OAT) sensitivity analysis and the Generalized Likelihood Uncertainty Estimation
(GLUE) was implemented for the two study watersheds. Results indicate that the suggested
range of some SWAT parameters, especially the ones that are used to determine the transport
capacity of channel network and initial concentration of nutrients in soils, required site-
specific adjustment. It was evident that uncertainties associated with sediment and nutrient
outputs of the model were too large, perhaps limiting its application for point estimates of
design quantities. However, the estimated effectiveness of BMPs sampled at different points
in the parameter space varied by less than 10% for all variables of interest. This suggested
that BMP effectiveness could be ascertained with good confidence using models, thus
making it suitable for use in watershed management plans such as the EPA's Total Maximum
Daily Load (TMDL) program. The potential impact of our analysis on utility of models and
model uncertainties in decision making process is discussed.
14
-------
3.2. Introduction
The analysis of uncertainty associated with the utility of simulation models is an important
consideration in the development of watershed management plans. Modeling uncertainty
should be rigorously addressed in development and application of models, especially when
stakeholders are affected by the decisions contingent upon model-supported analyses (NRC,
2001). Watershed models are commonly utilized to investigate rainfall-runoff generation, and
fate and transport of contaminants resulting from nonpoint source activities. Nonpoint source
activities are perceived to be the most important source of pollution in the United States (Ice,
2004). The evaluation of the success of best management practices (BMPs) in meeting their
original goals has also been facilitated by watershed models (Griffin, 1995; Edwards et al.,
1996; Mostaghimi et al., 1997; Saleh et al., 2000; Santhi et al., 2001; Kirsch et al., 2002;
Santhi et al., 2003; Arabi et al., 2006b). Uncertainty associated with absolute estimates of
design quantities tends to be very high because of data sparsity and model limitations
(Osidele et al., 2003; Benaman and Shoemaker, 2004). Thus, models are found to be more
useful when making relative comparisons rather than making absolute predictions. It may be
more meaningful to implement the uncertainty associated with effectiveness of BMPs in the
planning process.
The common modeling approach entails the (calibrate —>• validate —>• predict} process. The
thrust of the calibration procedure is to identify a set of model parameters by optimizing a
goodness-of-fit statistic between observed and predicted values such as the Nash-Sutcliff
coefficient of efficiency. The calibrated model is then used to examine the impact of various
management scenarios on the future behavior of the system. Such an analysis is subject to
identifiability, and non-uniqueness of the optimal (calibrated) parameter set (Beck, 1987), i.e.
there may be several sets of model parameters that fit the observed data equally (Beven and
Binely, 1992). Calibration of a simulation model for a given watershed will reduce, but not
totally remove, modeling uncertainties associated with both structure of the model and
parameter estimates. Even with the best model structure, parameter estimation contains
residual uncertainty (Beck, 1987) that propagates forward into model predictions and
evaluation of effectiveness of management practices.
15
-------
Although the literature is replete with sensitivity analysis and uncertainty analysis methods
(Spear and Hornberger, 1980; Beven and Binely, 1992; Spear et al., 1994; Saltelli et al.,
2000), implications of uncertainty associated with model predictions have not been widely
endorsed in the decision making process mainly as a result of large uncertainty estimates.
The magnitude of uncertainty itself is a key factor in its acceptance as the cost of
implementation of management actions such as the TMDL program may significantly
increase with larger uncertainty estimates (Dilks and Freedman, 2004). In a case study in the
Cannonsville Reservoir system watershed (1178 km2) located in upstate New York, Benaman
and Shoemaker (2004) concluded that even in the presence of observed data it was not
possible to reduce the uncertainty of absolute sediment predictions in their study to practical
values for the TMDL program. The argument, however, is that if the goal of a modeling
study is to examine the impact of management scenarios on water quality of a study area, it
may be neither practical nor necessary to incorporate large uncertainty of absolute
predictions in the decision making process. It would be perhaps more feasible (and more
desirable) to communicate and implement uncertainty of estimated effectiveness of
management scenarios rather than uncertainty of absolute predictions (Zhang and Yu, 2004).
Moreover, the importance of such a formulation would be particularly appreciated when the
reduction of a variable of concern (sediment, nutrients, etc) due to implementation of an
abatement action is less than estimated uncertainty of absolute predictions. In such cases,
evaluation of impact of management scenarios would not be inhibited by uncertainty of
model outputs.
The impact of modeling uncertainties on evaluation of management practices has not been
addressed sufficiently, as studies have generally focused on uncertainty of point predictions.
Specifically, a computational procedure than can be used to establish uncertainty bounds for
the estimated effectiveness of BMPs has not been developed to the best of our knowledge. In
this study, a Monte Carlo-based probabilistic approach is utilized (i) to develop a
computational procedure for analysis of uncertainty; (ii) to examine the effect of modeling
uncertainties on evaluation of long-term water quality impacts of BMPs using a distributed
watershed model, SWAT; and (iii) to provide a comparison between magnitudes of
uncertainties associated with absolute predictions versus effectiveness of BMPs. The analysis
16
-------
is demonstrated for two small watersheds in Indiana where water quality data were collected
and several structural BMPs were implemented.
3.3. Theoretical Considerations
3.3.1. Morris's One-at A Time sensitivity analysis
The OAT is a sensitivity analysis technique that falls under the category of screening
methods (Saltelli et al. 2000). In the OAT, each model run involves perturbation of only one
parameter in turn. This way, the variation of model output can be unambiguously attributed
to the perturbation of the corresponding factor. For each input parameter, local sensitivities
are computed at different points of the parameter space, and then the global (main) effect is
obtained by taking their average. The elementary effect of a small perturbation A of the /'th
component of the/>-dimensional parameter vector (a,) at a given point in the parameter space
a = (al,...,ai_l,ai,ai+l,...,ap)) is (Morris, 1991):
Eq31
where y(a) corresponds to model output. The results are quantitative, elementary, and
exclusive to the parameter a,. However, the elementary effect computed from Eq 3.1,
i.e. d(ai a) , is only a partial effect and depends on the values chosen for the other elements
of the parameter vector (o^). A finite distribution (F,) of elementary effects of parameter a,
is obtained by sampling at different points of the space, i.e. different choices of parameter set
a. The mean of the distributions is indicative of the overall influence of the parameter on the
output, while the variance demonstrates interactions with other parameters and nonlinearity
effects.
3.3.2. Generalized Likelihood Uncertainty Estimation (GLUE)
The GLUE methodology is based on recognition of the importance of the set of parameters to
produce the behavior of the system, not individual parameters. The acceptable model
realizations (i.e. behavior set) obtained from Monte Carlo simulations are given a likelihood
17
-------
weight according to observed data and a likelihood function. Several choices for an
appropriate likelihood function can be obtained from Beven and Freer (2001). A likelihood
measure based on Nash-Sutcliff efficiency criterion with shaping factor TV is defined as (Freer
etal., 1996):
Eq3.2
cr
where L(a\y) is the likelihood of parameter set (a), given the observed data (y). The quantities
cr? and &1 refer to the error variance between model simulations and observed data, and the
variance of the observed data, respectively. For N=l, L in Eq 3.2 is the well-known Nash-
Sutcliff efficiency coefficient that is often used for calibration of hydrologic and water
quality models. A negative L value indicates that the corresponding model output is
dissimilar to the behavior of the system under study, and the likelihood of such a simulation
in mimicking the system behavior is zero. Therefore, the likelihood measure is rewritten as:
Eq3.3
For each simulation from a random parameter set, a likelihood weight is obtained from Eq
3.3. Then, these weights are rescaled by dividing each of them by their total sum. The
rescaled likelihood weights are used to construct a cumulative distribution for the output of
interest, which can be used for estimation of uncertainty bounds associated with the output
by computing its quantiles. The GLUE method is subjective not only to the choice of
likelihood function, but also to the cutoff criterion used for classification of Monte-Carlo
simulations to behavior and non-behavior sets.
3.4. Methodology
The computational framework utilized in this study is comprised of (1) a process-based
watershed model (SWAT; Soil and Water Assessment Tool, Arnold et al., 1998) employed
for simulating the fate and transport of sediment and nutrients in the study watersheds under
two scenarios- before and after implementation of BMPs, (2) a BMP representation method
18
-------
(adapted from Arabi et al., 2004) for incorporation of water quality impacts of BMPs, and (3)
an uncertainty analysis methodology based on two sampling based methods: One-factor-At-
a-Time (OAT) sensitivity analysis (Saltelli et al., 2000), and Generalized Likelihood
Uncertainty Estimation (GLUE; Beven and Binely, 1992). A computational method was
developed to reduce the "suggested" range of input parameters to site-specific "adjusted"
ranges to be used in the GLUE analysis. The methodology was tested on two different
watersheds located in northeast Indiana, one with significantly larger number of installed
BMPs. This confirmed the versatility of the computational analysis to quantify modeling
uncertainties associated with estimated effectiveness of BMPs in study watersheds under
influence of a large number as well as a small number of BMPs.
3.4.1. OAT sensitivity analysis
Identifying the most sensitive input parameters is not an essential component of the
framework for analysis of uncertainty that is proposed in this study. The results of such an
analysis, however, provide useful information pertaining to model parameters with large
uncertainties. The uncertainty associated with parameters that are determined based on
topographic, land use, soil, and management attributes could perhaps be reduced by using
better spatial resolution for these attributes. Conversely, uncertainty of parameters that are
not determined from landscape attributes and are usually estimated through calibration
procedure may be reduced only through a systematic site-specific range adjustment process.
In this study, the following procedure was followed for the OAT sensitivity analysis. First, in
the absence of any prior information with regard to the distribution of input parameters, a
uniform distribution, with minima and maxima specified in Table 2.2, was considered for all
parameters. The same approach was used by Beven and Freer (2001) and Benaman and
Shoemaker (2004). Then, OAT sensitivity indices were computed. While a local sensitivity
index (d) was computed by applying Eq 3.1 for each parameter, global sensitivities were
determined by taking the average of these local sensitivities at 50 random points sampled
across the entire parameter space. Finally, a rescaled sensitivity index (dn~) was determined by
dividing the global OAT sensitivity indices by their total sum. The dn values range from
[0,1], with values closer to 1 indicating a more sensitive parameter.
19
-------
3.4.2. Range adjustment
Adjustment of "suggested" ranges of sensitive input parameters was found to be an essential
step prior to quantifying uncertainty of model outputs. Use of unrealistically large parameter
ranges may result in a large number of simulations with negative likelihoods (L in Eq 3.2)
that renders the GLUE analysis inefficient and prohibitively expensive. Particularly, ranges
of calibration parameters that are not determined from landscape attributes should be
investigated. For example, consider the model parameters that are used to determine the
transport capacity of the channel network (SPCON, SPEXP, and PRF in Table 2.2). These
parameters are not identified based on topographic and/or soil characteristics of the channel
segments, but typically determined based on similar research in the study area or through a
calibration procedure.
The range adjustment process for important model parameters entailed the following steps:
1. Establish base-case values for input parameters: the base-case values for the input
parameters in Table 2.2 could be selected by one the following approaches: (i) insights
gained from the OAT sensitivity analysis; (ii) model calibration procedure; (iii) a
suggested value obtained from literature, previous studies in the study area, or prior
experience of the analyst. In this analysis, base-case values were selected from a manual
calibration (Arabi et al., 2004).
2. Perform interval-spaced simulations: the range of each parameter of concern was divided
into 50 equal intervals. A SWAT model simulation was performed for a random
realization of the parameter from each interval while other parameters were kept at their
base-case values. The simulation timeline covered the 1974-1978 period when hydrologic
and water quality data were collected at the outlets of the study watersheds.
3. Adjust ranges of important parameters: A method based on a goodness-of-fit measure
(Nash-Sutcliff coefficient) was developed to reduce the suggested ranges of important
parameters to narrower adjusted ranges for the study watersheds. The acceptable ranges
of parameters were determined by plotting the Nash-Sutcliff efficiency coefficients (E^.s,
Nash and Sutcliff, 1970) computed for interval-spaced simulations. The E^.s is defined
as:
20
-------
where Oi and /J are observed and predicted monthly output variables, respectively, while
O is the average of monthly observed values. EN-S ranges from -oo to 1, with higher
values indicating a better 1:1 fit between the observed and simulated values. Model
simulations with negative EN-S were considered "unacceptable" (Santhi et al., 2001), and
subsequently, the range of the parameter was reduced to those that yielded non-negative
EN_s values.
3.4.3. Uncertainty computations
Adjusted parameter ranges from previous step were used to probabilistically determine the
effectiveness of BMP implementation. In the present uncertainty analysis, 5000 model
simulations were performed for the 1974-1978 period when water quality data were
collected. All model parameters in Table 2.2 were included in the analysis. For those
parameters that required range reduction, adjusted ranges were selected from the range
adjustment process. For each realization of the parameter space, the uncertainty analysis
entailed the following steps:
1. Select parameter values from their specified ranges in a random fashion;
2. Perform model simulation with the selected values to compute model outputs for the
scenario without BMPs;
3. Compute likelihood of the simulation for monthly output variables by applying Eq 3.3;
4. Represent BMPs by changing appropriate model parameters from Table 2.3;
5. Compute model outputs for the scenario with BMPs in place;
6. Using GLUE-likelihood weights computed at step (3), generate a cumulative density
function for outputs of interest for both model simulations with and without BMPs.
The main assumption in the numerical procedure described above is that the same cumulative
likelihood can be used for the scenarios with BMPs and without BMPs. Water quality data
utilized in this study were representative of the state of the watersheds before implementation
21
-------
of BMPs. Thus, the likelihood of a set of model parameters to mimic the behavior of the
system was computed for the scenario without BMPs. When adequate hydrologic and water
quality data are available for both scenarios with BMPs and without BMPs, the GLUE-
likelihood functions should be computed separately.
It should be noted that likelihood functions associated with each set of model parameters in
the GLUE analysis are likely to be approximately the same before and after implementation
of BMPs, if other watershed characteristics remain the same. Representation of BMPs entails
only modification of some model parameters for the fields and/or streams where BMPs are
implemented. BMPs are implemented in a relatively small area within the watershed. For
example, BMPs in Figure 2.1 cover less than 5% of total upland fields and the channel
network in the Dreisbach and Smith Fry watersheds. Thus, in the absence of water quality
data for computation of GLUE-likelihoods for the scenario with or without BMPs, it is
suggested that the same likelihood be used for both scenarios.
3.5. Results
3.5.1. Range adjustment results
The results of the range adjustment analysis indicated that none of hydrology-related
parameters of the SWAT model needed range reduction. On the other hand, three sediment-
related parameters, one total P-related parameter, and one total N parameter required site
specific adjustments. Table 3.1 provides a listing of top five sensitive SWAT parameters in a
descending order for monthly streamflow, sediment, total P, and total N output variables. The
term "min (Ejv-s)" in the table refers to the minimum of Nash-Sutcliff coefficients computed
for 50 interval-spaced model simulations corresponding to each parameter. Shaded values
indicate parameters with negative "min (Ejv-s)" value for which the range adjustment
procedure was applied.
22
-------
Table 3.1 Top five sensitive SWAT parameters for streamflow,
sediment, total P, and total N computations: dn is the OAT sensitivity
index, "min (EN_s)" refers to minimum of Nash-Sutcliff coefficient
computed for each 50 simulations in the OAT analysis. Description
of parameters can be found from Table 2.2.
Variable
Streamflow
Sediment
Total P
Total N
SWAT Symbol
SOL_AWC
CN
GW_REVAP
ESCO
CH_KI
SPCON
PRF
SLOPE
CH_N2
CN
SLOPE
ORGP_AG
LABP_AG
CN
USLE_K
SLOPE
ORGN_AG
USLE_K
CN
SOL_AWC
dn
0.30
0.20
0.17
0.09
0.09
0.30
0.12
0.09
0.06
0.06
0.43
0.25
0.07
0.06
0.05
0.51
0.11
0.06
0.07
0.06
min
Dreisbach
0.19
0.38
0.27
0.61
0.40
-30.59
-0.31
-2.07
-0.27
0.60
-110.02
-1.69
0.03
0.09
0.12
-132.94
0.03
0.12
0.18
0.14
(EN.s)
Smith Fry
0.13
0.13
0.36
0.56
0.60
-4.9
-0.85
0.11
-0.2
0.01
-10.66
-0.19
0.4
0.12
0.05
-11.77
-0.5
0.05
0.05
0.34
Figure 3.1 illustrates one dimensional response surface of model outputs at the outlets of the
study watersheds to parameters with reduced ranges. The spider-plots (dashed lines) have
been generated by varying one parameter at a time, while other parameters were kept
constant at their base-case value. Each panel represents 50 model simulations, where the
model output is shown on the left y-axis and the right y-axis displays the Nash-Sutcliff
coefficient in Eq 3.4. The x-axis is the normalized value of each parameter a., determined
from its absolute valueai, and its respective upper (t/) and lower (Z.) bounds summarized
in Table 2.2:
23
-------
Dreishach Watershed
0.2 0.4 0.6 0.8
SPCON
I
-£± o
C f-
E ^
5 <£-
0.04
0,03
0.02
0.01
o
0.04
0.03
0.02
0.01
0
0
0.2 0.4 0.6 OJ
pfff
0
02 0.4 0.6 0.8
CH N2
ORGP,
1
1
1
075
0.5
025
0
1
OS
o ,
-0.5
-1
• Average Monthly Output Variable
"WS
036-
Smith Fry Watershed
!-. 0.27-
£1
II 0.18'
E 3
I e 0.09 •
0,12
1 0.09
>-£
1 I 0,06
E w
TjJ=
S - 0,03
0
0,12
S? 009
>f
| | 0,06
E re
•a€
0.2 0.4 0.8 0.8
SPCO/V
0.2 0.4 0.6 0,8
PRF
0.2 0.4 0.6 0.8
CH N2
0.76
0.2 0.4 0.6 0.8
QRGf'
fi 4,5
! c
• o
ic
;f 15
0
0.2 0.4 0.6 0.8
-05
-2
-3.5
1
0.66
0.33
0
-0.33
1
07
0.4 u
0.1
-0.2
1
1
0.5
W3
0 uj*
-0,5
1
Figure 3.1 Spider plots for the most sensitive input factors related to water quality
computations in SWAT with adjusted ranges in Table 3.2. Each panel represents 50 model
realizations. Subscript "AG" refers to agricultural land use.
24
-------
-
a, = •
Eq3.5
Equation 3.5 can be used to back calculate the absolute parameter values corresponding to
the normalized values shown in Figure 3.1. A negative "min (EN_s)" value indicated the
necessity for range reduction. For each panel, a horizontal line was drawn at EN-S (right y-
axis) equal to zero. The range of the corresponding parameter vector was reduced to the
portion that lies above this line. The top panel to the left in Figure 3.1 demonstrates how
parameter ranges were adjusted. The new parameter ranges for the study watersheds are
presented in Table 3.2.
Table 3.2 Adjusted parameter ranges for the study watersheds: both suggested ranges and
adjusted ranges are based on absolute parameter values. Description of the parameters can be
found in Table 2.2.
Parameter Units Limiting Variable Suggested Range
Adjusted Range
Dreisbach Smith Fry
SPCON
PRF
CH_N2
ORGP AG
ORGN AG
mg/kg
mg/kg
sediment
sediment
sediment
total P
total N
0-0.001 0-0.0002 0-0.0005
0-2 0.2-2 0.2-2
0.008-0.3 0.008-0.1 0.008-0.1
0-1000 0-500 0-950
0-10000 0-10000 2000-10000
Sediment outputs were most sensitive to model parameters that are used to estimate the
transport capacity of the channel network, usually determined through calibration procedure
and not based on land and/soil/management characteristics. The slope of upland areas
represented by parameter SLOPE was the most sensitive parameter for total P and total N
simulations. However, the range of this parameter was not reduced, because it is directly
estimated from Digital Elevation Model (DEM). In the absence of field measurements, the
initial concentration of phosphorus and nitrogen in soils were reduced to the values reported
in Table 3.2.
Several issues should be considered in interpretation of the results from range adjustment
analysis. First, the threshold value of Nash-Sutcliff coefficient that is used to reduce
25
-------
parameter ranges significantly impacts the method. Obviously, a lower EN-S value would
result in acceptance of wider parameter ranges. In this study, this threshold was set at zero
because the likelihood of model simulations in GLUE uncertainty analysis with negative EN_s
values is assumed to be zero (see Eq 3.3). Second, the number of intervals to be used in the
range adjustment process should be selected such that the plot of the goodness-of-fit measure
(EN-S) over the entire parameter range (e.g. Figure 3.1) forms a relatively smooth curve. Our
experience from this study indicated that 20-50 intervals should be adequate. The results of
the OAT sensitivity analysis and range adjustment process are site specific, and may be
different if the watershed conditions are significantly different than those discussed here. We
note that the results of the sensitivity analysis are in general agreement with previous studies
such as Lenhart et al. (2002), Heuvelmans et al. (2004), Muleta and Nicklow (2005), and
White and Chaubey (2005).
3.5.2. Uncertainty analysis results
The computational procedure described previously was performed to quantify the uncertainty
associated with both model simulations and estimated effectiveness of BMPs. Figures 3.2
and 3.3 depict the cumulative likelihood distribution of average monthly values for various
output constituents at the outlet of the Dreisbach and Smith Fry watersheds, respectively. A
likelihood weight from Eq 3.3 was assigned to each of 5000 model realizations resulting
from Monte Carlo simulations. Each panel in the figures contains two cumulative probability
distributions for the two scenarios. The dashed lines correspond to results from the GLUE
method for the scenario without BMPs (scenario A), while the solid lines demonstrate the
results for the scenario with BMPs represented in the model (scenario B). Two major trends
are evident. First, the difference between GLUE-likelihoods for scenarios A and B was
marginal for streamflow, but substantial for sediment, total P, and total N computations.
Second, the trends were quite different in the two study watersheds.
26
-------
0.81-
0.61
o I
f 0.4J
§ 0.21-
d i
BMP* Not Represented
• BMPs Represented
0.025
0.05
0.075
0.1
Monthly Gtreamflow (m /s)
BMPs No! Represented
BMPs Represented
0.09
0.12
Yi*Md {
BMPs Not Represented
BMPs Represented
0.05 0.1 0,15
Average Monthly Tola! P Yield (kg/ha/month)
0.2
BMPs Not Represented
BMPs Represented
0 0.75 1.5 2.25
Average Monthly Total N Yield Skg/ha'moniti)
Figure 3.2 Cumulative GLUE-likelihood for variables simulated at the outlet of the
Dreisbach watershed, Indiana. Each panel represents 5000 model realizations. The top panel
in the right demonstrates how various percentiles were determined. Variable ya is the
expected value for scenario A (BMPs not represented), and ybis the expected value for
scenario B (BMPs represented).
27
-------
0,8
u 0.6 '
3 i
d I
f 0.4 [
0.2 r
BMPs Not Repre»ented
• BMPs Represented
0.03 0.06
AvBra^jFi Mnntfily Str
0.09
m*/
0.12
OMPs Not Represented
BMPs Represented
0 0.03 0.06 0.09
Average Monthly Sediment Yield (t/hai'montrt)
0.12
BMPs Not Represented
BMPs Represented
Average Monthly Total P Yield (kg/ha/morsth)
BMPs Not Represented
BMPs Represented
0248
Average Monthly Total N Yield (kgflia/montti)
Figure 3.3 Cumulative GLUE-likelihood for variables simulated at the outlet of the Smith
Fry watershed, Indiana. Each panel represents 5000 model realizations.
28
-------
Following Figures 3.2 and 3.3, the results are summarized in Table 3.3 for further discussion.
Lower and upper bounds in Table 3.3 refer to 5th and 95th percentiles of the cumulative
likelihoods, respectively, while the median represents the 50th percentile. The top panel on
the right in Figure 3.2 demonstrates how various percentiles were determined. For scenario
A, comparison of median values to the average of monthly observed data, also presented in
Table, indicate that averages of monthly observed values for streamflow were within ±15%
of the median in both study watersheds, and fell well within the uncertainty ranges from
GLUE. Sediment yield measurements at the outlet of the Dreisbach and Smith Fry
watersheds were within ±20% of the median values, also well covered by the uncertainty
bounds. The median values for total P appeared to adequately match the average of monthly
observations as they underpredicted by only 27% and 20% in Dreisbach and Smith Fry,
respectively. Moreover, the uncertainty bounds covered the data. The median for total N
determined for Dreisbach watershed differed from the average of monthly observed values
by only 7%, indicating satisfactory results from GLUE. Total N at the outlet of Smith Fry
was the only case where the uncertainty bounds from GLUE did not contain the observed
value. The corresponding median underpredicted total N yield by nearly 50%, perhaps
because of structural uncertainties involved in simulations. In particular, a lack of proper
linkage between fluvial channel processes and in-stream nutrient processes in SWAT could
be responsible for these results.
Table 3.3 Summary of the results from analysis of uncertainty: variable y refers to 50l
percentile (median) of output variables for 5000 Monte Carlo simulations of the SWAT
model, f is the 50th percentile (median) of estimated reduction of output variables due to
implementation of BMPs computed by 3.3 for 5000 Monte Carlo simulations of the SWAT
model, and "range" refers to the 5th and 95th percentile of the Monte Carlo simulations.
Watershed Variable
/- Streamflow
o
,g Sediment
'| Total P
0 Total N
^ Streamflow
^ Sediment
•g Total P
00 Total N
Scenario A
Units Observed Without BMPs
Value
ya range
nrVs
t/ha/month
kg/ha/month
kg/ha/month
nrVs
t/ha/month
kg/ha/month
kg/ha/month
0.039
0.031
0.095
1.228
0.054
0.093
0.336
6.112
0.041
0.035
0.069
1.145
0.06
0.043
0.267
3.026
[0.026,0
062]
[0.015,0.083]
[0.030,0
[0.487,2
[0.045,0
[0.021,0
[0.115,0
[1.552,5
125]
385]
086]
107]
551]
838]
Scenario B
With BMPs
0.037
0.015
0.046
0.788
0.055
0.037
0.243
2.847
[0.023,0
[0.007,0
[0.015,0
[0.360,1
[0.042,0
[0.018,0
[0.098,0
[1.450,5
056]
033]
086]
502]
077]
090]
502]
560]
Reduction (%)
ft range
9.8
57.1
33.3
31.2
8.3
14
9
5.9
[9.70
[53.3
[31.2
[26.1
[6.70
[12.8
[8.90
[4.70
,11.5]
,60.2]
.4,50]
,37.0]
,10.5]
,15.9]
,14.8]
,6.60]
29
-------
These values may now be compared to corresponding quantities for scenario B listed in
Table 3.3. Comparison of the expected values for scenarios A and B revealed the estimated
effectiveness of BMPs, shown as a percent reduction. The percent reduction of each
constituent as a result of implementation of BMPs was determined as:
~ lOO Eq3.6
yt,a
where /' is the constituent of interest, i.e., streamflow, sediment yield, total P load, or total N
load, rt is percent reduction of constituent /', y^a is the median of constituent /' for scenario A,
andyib is the median of constituent /' for scenario B. The results indicated that
implementation of BMPs as shown in Figure 2.1 would not reduce streamflows significantly.
This did not come as a surprise because the BMPs implemented in the watersheds were
mostly sediment control BMPs. Implementation of the 26 BMPs in the Dreisbach watershed
would lower median sediment yield by nearly 57%, from 0.035 t/ha/month to 0.015
t/ha/month. The estimated reduction rates for total P and total N at the Dreisbach outlet were
33% and 31%, respectively. The estimated reduction rates at the outlet of Smith Fry were
nearly 14%, 9%, and 6% for sediment yield, and total P and total N loads (Table 3.3). These
rates suggested that implementation of the BMPs in Smith Fry was almost four times less
effective than the ones in Dreisbach, which was anticipated because of smaller number of
BMPs in the former.
The uncertainty bounds associated with absolute predictions of the SWAT model were much
larger than the ones corresponding to the estimated effectiveness of BMPs. For example, the
uncertainty bound for monthly sediment yield at the outlet of Dreisbach for the simulation
period under scenario A (scenario with no BMP represented) was determined to be
[0.015,0.083] with a median value of 0.035 t/ha/month (Table 3.3). Incorporation of this
large uncertainty in decision making and management may be extremely costly and not
feasible. However, as is evident in Table 3.3, the estimated 5th and 95th uncertainty bounds
for estimated effectiveness of BMPs in the Dreisbach watershed in reducing sediment yield
at its outlet were [53%,60%] with a median (50th percentile) of 57%. Likewise, the estimated
30
-------
effectiveness of BMPs in Smith Fry was estimated to range between [13%, 16%] with a
median of 14%. Similar results were observed for streamflow, total P, and total N
computations for the study watersheds.
3.6. Discussion
The conceptual simplicity is an attractive feature of the computational analysis developed in
this study. The utilized likelihood measure, Nash-Sutcliff coefficient, is widely used by
modelers for calibration and validation of watershed models. Also, the OAT-GLUE
methodology will, more than likely, produce various sets of model parameters that
adequately satisfy calibration criteria that are usually set based on Nash-Sutcliff coefficient.
This helps modelers avoid the cumbersome practice of manual calibration. That the
uncertainty bounds of various constituents encompassed the corresponding observed values
(except for only one case for total N in the Smith Fry watershed) strongly supports this
suggestion. It appears that for the Dreisbach and Smith Fry watersheds, sediment and nutrient
outputs of the SWAT model bear more uncertainty than streamflow simulations.
Consolidation of results in Tables 3.1 and 3.3 from the three-step procedure indicate that
reducing the uncertainty associated with absolute sediment and nutrient outputs of the SWAT
model to practical ranges may not be feasible. Sediment outputs of the model were found to
be most sensitive to transport capacity of channel network. Some of the parameters that are
used to determine the transport capacity of channel segments (SPCON, SPEXP, and PRF in
Table 2.2) cannot be measured in the field and are usually calibrated from a broad
"suggested" range. With availability of more data, the three-step procedure used in this study
may result in more narrow uncertainty bounds for absolute predictions. Nevertheless, it may
still not be adequate for reducing the uncertainty associated with these absolute model
predictions to small enough ranges to be useful for practical management decisions.
However, comparison of sediment and nutrient outputs from model simulations with and
without representation of BMPs would not suffer from such limitations as demonstrated here.
It is evident that incorporation of modeling uncertainty in informed decision making is more
practical through communicating the uncertainty associated with evaluation of effectiveness
of management actions.
-------
The foregoing results indicate that the computational framework developed in this study
could be endorsed by watershed management programs such as the TMDL as advocated in
NRC (2001). Once the level of protection to be provided by the Margin of Safety (MOS) is
specified by decision makers, the probabilistic framework presented in this study could be
applied to verify the success of a particular management action in achieving its designated
goals. Implementation of the estimated MOS is more feasible in this context, as opposed to
the MOS computed from uncertainty analysis of absolute model predictions where the
additional cost of implementation may render its consideration impractical. Finally, the
methodology can be easily incorporated in watershed models such as SWAT that are
commonly used for TMDL development. Land use in the study watershed has changed since
data were collected. Thus, application of the results presented herein is limited to
demonstration of the methodology.
3.7. Conclusions
A computational framework was developed in which investigation of uncertainty provides
complementary quantitative and qualitative information to support management and decision
making. The analysis focused specifically on two issues. First, ranges of some model
parameters may require site-specific adjustments with available data. This was particularly
evident for parameters that are not determined from landscape characteristics. Second, the
uncertainty associated with estimated effectiveness of BMPs is substantially smaller than the
uncertainty associated with absolute predictions. The computational procedure for analysis of
uncertainty was performed for two Indiana watersheds with relatively similar spatial scale,
and landscape characteristics, but different number of installed BMPs. It was demonstrated
that the computational framework is capable of quantifying uncertainty of effectiveness of
implemented BMPs in both watersheds. Future work will focus on investigation of the utility
of the developed methodology in the TMDL process. Especially, the means for coping with
the Margin of Safety (MOS) in the TMDL formula should be explored.
32
-------
SECTION 4. SENSITIVITY ANALYSIS OF SEDIMENT AND NUTRIENT PROCESSES WITH A
WATERSHED MODEL
4.1. Abstract
This section presents a computational analysis for evaluating critical nonpoint source
sediment and nutrient processes and management actions at the watershed scale. In the
analysis, model parameters that bear key uncertainties are presumed to reflect the importance
of natural processes and/or management actions that they represent. The Regionalized
Sensitivity Analysis (RSA) and the Tree-Structured Density Estimation (TSDE) procedures
were integrated with the Soil and Water Assessment Tool (SWAT) to investigate correlation
structure in the parameter space while accounting for multiple objectives. The computational
analysis was applied to the Dreisbach and Smith Fry watersheds in Indiana in the Midwestern
portion of the United States. Results showed that incorporation of parameter interactions is
essential to obtaining conclusive information about critical system processes and
management actions. Interactions between surface runoff volume and within-channel
processes were critical to describe transport of sediments in the study watershed. Key
management actions for nutrient control were found to be fertilizer application and upland
farming practices such as parallel terraces. The sensitivity analysis reported herein could be
used to derive a list of key nonpoint source best management practices for development of
watershed management plans.
4.2. Introduction
Identification of natural processes and management actions that control nonpoint source
(NFS) pollution is essential for development of watershed management plans. The current
handbook of the Environmental Protection Agency for developing watershed plans (EPA,
2005a) recommends implementation of best management practices (BMPs) in portions of the
watershed that are believed to contribute intensively to NPS pollution. Locating critical areas,
however, is complicated because contaminants are carried along with the flow, and water
33
-------
movement over a watershed tends to be fairly dynamic, behaving in a nonlinear fashion.
Potential interactions among hydrologic, fluvial, and nutrient processes and NFS control
measures should be identified and included in designing watershed management scenarios.
Such information could be obtained by integration of computational techniques with
hydrologic/water quality models.
Previous studies have utilized optimization methods to design near optimal combination of
BMPs (Srivastava et al., 2002; Veith et al., 2004). These methods search for the optimal
location of selected management practices for NFS pollution control in a watershed system.
As the number of selected BMPs increases, the number of model evaluations for identifying
the optimal solution increases exponentially. Thus, these optimization studies have focused
on spatial allocation of only a few BMPs. Nearly a hundred BMPs for NFS pollution control
are recommended by the USDA Natural Resources Conservation Service (USDA NRCS,
2006). Evaluation of impacts of all these practices, even within small watersheds, is
infeasible. In this section, we present a novel approach to systematically abridge the list of
key BMPs for sediment and nutrient control at the watershed scale. The methodology is
based on the hypothesis that analysis of uncertainty of model simulations could underscore
critical processes and key management actions for NFS control.
Inverse modeling approaches based on sensitivity analysis (SA) have been used in the past to
obtain information about important system processes in a variety of disciplines. Osidele et al.
(2003) utilized a sensitivity analysis to investigate the importance of sediment and nutrient
processes in a section of the Chattahoochee River (nearly 115 river miles) in Atlanta. The
univariate regionalized sensitivity analysis (RSA; Spear et al., 1980) was used in conjunction
with the multivariate tree structured density estimation (TSDE; Spear et al., 1994) to account
for the correlation structure in the parameter space. The RSA was also utilized by Zheng and
Keller (2006) for sensitivity analysis of the Watershed Analysis Risk Management
Framework (WARMF) model and its management implications.
The one-at-a time sensitivity analysis (Pitman, 1994), factorial design, and the Fourier
amplitude sensitivity test (FAST) (Saltelli et al., 2000) are some other SA methods. These
SA methods determine parameter sensitivities for a single model response (e.g. average
34
-------
annual sediment yield). Bastidas (1998) developed a multiobjective generalized sensitivity
analysis (MOGSA) based on RSA to consider the influence of multiple criteria on parameter
sensitivities. Previous research studies based on MOGSA (Bastidas et al., 1999; Meixner et
al., 1999; Liu et al., 2004, Demarty et al., 2005), while incorporating multiple criteria, were
limited in revealing the multivariate correlation structure in the parameter space. Moreover,
these studies have been applied only to models up to medium level of complexity (Zheng and
Keller, 2006).
The work of Osidele et al. (2003) and van Griensven et al. (2006) are examples of studies
that applied multivariate approaches for sensitivity analysis of sediment and nutrient
processes. However, they evaluated model sensitivities based on a single output or its
trajectory, and did not explicitly address the interactions between flow, sediment, and
nutrient criteria. To our knowledge, both multivariate interactions and multiple criteria have
not been incorporated simultaneously for performing a sensitivity analysis with a complex
watershed model such as SWAT (Arnold et al., 1998). Also, implications of sensitivity of
NFS sediment and nutrient processes/BMPs in holistic watershed management have not been
discussed. These shortcomings are addressed in the present work.
In this section, we investigate the utility of computational approaches for identifying critical
NFS processes and key BMPs that are likely to control fate and transport of sediments and
nutrients in watershed systems. To this end, the RSA and TSDE procedures are linked with a
comprehensive watershed model (SWAT) to (i) develop a new sensitivity analysis
framework that can handle interactions in the parameter space as well as multiple trajectories
of model outputs, and (ii) demonstrate the application of the computational analysis to
highlight critical sediment and nutrient processes in an agricultural watershed in Indiana.
With simultaneous incorporation of multicriteria and multivariate approaches, the developed
framework provides a novel capability for dealing with important issues for holistic
watershed management. Critical NFS sediment and nutrient processes, key agricultural
BMPs, and their potential interactions can be evaluated. Watershed management programs
such as the Total Maximum Daily Load (TMDL) and the USDA's Environmental Quality
Incentive Program (EQIP) would significantly benefit from such an analysis. Another feature
of the present work is that real system response data are utilized in the analysis.
35
-------
4.3. Theoretical Considerations
4.3.1. Regionalized Sensitivity Analysis (RSA)
The Regionalized Sensitivity Analysis (RSA; Spear and Hornberger, 1980) aims at
identification of critical uncertainties in order to evaluate the relative importance of
individual parameters that exert the most influence on system behavior. The procedure
utilizes a uniform sampling of the parameter space and involves (i) a qualitative definition of
the behavior of the system under study (criterion function C), and (ii) a binary classification
of the parameter space (S) into good (behavior^) or bad (nonbehavior 5) regions. The
strength of the method lies in the classification scheme, which facilitates the application of
multivariate statistical methods to explore the level of significance that the posterior
probability distribution function of each element of the parameter vectors in the behavior
region/m(dr B) deviates from the one in the nonbehavior region/n(d? 5). A Kolmogorov-
Smirnov two sample test is performed to test the hypothesis that for a given parameter
ak e d, fm (ak B) separates from /„ (ak B):
H0: fm(ak B} = fn(ak\B)
H,: fm(ak B)*fn(ak B) Eq4.1
Kolmogorov-Smirnov statistic: dmn= sup Fm(ak B)-Fn(ak B}\
X
where Fm(ak B) and Fn(ak B) are, respectively, the cumulative density functions
corresponding to fm(ak B) and fn(ak B) for m behavior and n nonbehavior model
simulations. The sup notation refers to the largest vertical separation between Fm(ak B)
x
andFn(ak B). The dmnstatistic can be used to determine the relative importance of the
uncertainty associated with each element of the parameter vector, with higher values
indicating higher influence on model outputs.
4.3.2. Tree Structured Density Estimation (TSDE)
Spear et al. (1994) recognized that despite the conceptual simplicity of the RSA procedure,
its applicability may be limited as a result of different correlation structures among
36
-------
parameters. The tree structured density estimation (TSDE) methodology was developed to
obtain information relevant to the interactions between model parameters in the complex
non-uniform behavior region. The TSDE procedure is as follows. The m behavior parameter
vectors are treated as independent samples from an unknown probability distribution
function /(a 1 5). To construct an adequate approximation of this unknown distribution, the
behavior parameter space (SW) is partitioned into q sub-spaces (SB = \ \q SBi\ A local
>-«'!=l
estimate off for the ith sub-space (83,1) with volume (Vj) can be defined as (Spear et al., 1994):
-x() Eq4.2
m V
where mi is the number of behavior parameter vectors in SB,I. With/? reflecting the number of
individual model parameters, the [m x p] parameter hypercube can be split into q sub-spaces
in q x p ways. The search algorithm for finding the optimal split involves minimizing a loss
function (L) that mimics deviation of f(ak \ B) from f(ak B) and can be defined as (Spear et
al., 1994):
Ecl4-3
If the parameters were uniformly distributed in the parameter space SB (i.e. no split was
required), the first estimated density function (/0) and the first loss function (Z0) would be
equal to l/F and-l/F, respectively, where V reflects the volume of the behavior space (B).
The splitting process is performed successively, and begins with splitting the parameter space
into two sub-spaces. While the loss function corresponding to the first split L\ is always less
than LQ, LQ- LI is a measure of accuracy of approximation of the density function. Therefore,
maximizing LQ- L\ is synonymous with minimizing the errors associated with density
estimation. The parameter space is split on the axis of the parameter that produces the largest
increase in the accuracy criterion (L0- LI). Likewise, in the second recursion, each of the two
sub-spaces constructed in the first step is split into two new sub-spaces. This procedure is
37
-------
repeated until either the accuracy of density estimate does not increase significantly or the
density of each sub-space is less than some critical value.
The TSDE procedure finally yields a tree structure that reveals the multivariate correlations
among model parameters. The top (origin) node in the tree represents the original sample
space with normalized relative density, defined as the volume of each sub-space to the
volume of the original pace, equal to 1. The density of end (terminal) nodes indicates the
relative importance of the corresponding intermediate parameters and their interactions for
matching the behavior of the system under study. Beginning from the origin node and ending
with a terminal node, each branch graphically depicts the interactions between model
parameters.
4.4. Methodology
The computational analysis presented in this section was aimed at identification of
hydrologic and water quality processes that are likely to control transport of sediments and
nutrients. Model parameters served as surrogates for processes that they represent in the
model's mathematical structure. The resulting ranking and classification of the parameters
based on the importance of their uncertainties indicated the relative importance of the system
processes they represent. Inferences relevant to key management actions for sediment and
nutrient control were drawn.
The computational procedure shown in Figure 4.1 was implemented as follows. First, the
natural process and/or management action represented by each input parameter of the SWAT
model was identified. Model parameters and their suggested ranges are presented in Table
2.2. For example, CH_N2 represents the impact of roughness of the channel network on
sediment and nutrient transport. Likewise, USLE P indicates the importance of upland
farming practices such as parallel terraces (Renard et al., 1997). Next, 5000 parameter
vectors were randomly generated with a Latin Hypercube Sampling (LHS; McKay et al.,
1979) strategy. The SWAT model was used to simulate monthly flow, sediment and nutrient
outputs corresponding to each parameter vector. Model inputs and outputs were integrated
with the RSA and the TSDE methods to investigate significance of input factors.
38
-------
Input Factors
Natural / Fluvial
Processes
Management Actions
Random Sampling
Latin Hypercube
Sampling
(LHS)
Hydrologic / Water
Quality Model
Soil and Water
Assessment Tool
(SWAT)
Sensitivity Analysis
Behavior
Definition
0
Regionalized
Sensitivity Analysis
(RSA)
Tree-Structured
Density Estimation
(TSDE)
Key Input Factors
Critical Natural /
Fluvial Processes
Key Management
Actions
Figure 4.1 Computational framework for sensitivity analysis of sediment and nutrient
processes.
4.4.1. Behavior Definition
The definition of the behavior of a given system refers to a combination of thresholds,
extremes, and time scales derived from available or proposed information about the system
conditions. This combination defines a corridor through which the model output trajectory
must pass in order to qualify as a behavior simulation (Beck, 1987). Depending on the goals
of the study, any of the following can define the system behavior:
analysis of observed data,
speculations with regard to future state of the system,
• desired state of the system based on regulatory standards, and
target values for TMDL implementation.
The definition of the behavior of a given watershed system depends on the goal of the
problem at hand (Bastidas et al., 1999). In the present work, model simulations were
compared to the observed data and classified as behavior or non-behavior based on the Nash-
Sutcliff efficiency coefficient EN-S. Available data for the watersheds can be obtained from
Arabi et al. (2004). For streamflow, EN-S associated with streamflow output should be greater
39
-------
than or equal to zero in order to classify the simulation as behavior. For sediment yield, EN-S
associated with both streamflow and sediment output should be greater than or equal to zero
for a behavior simulation. A simulation was deemed behavior for the total N constituent only
if the EN-S values associated with all three streamflow, sediment, and total N outputs were
greater than or equal to zero. The Nash-Sutcliffe coefficient (Nash and Sutcliffe, 1970) that is
commonly used for calibration of hydrologic models can be defined as:
= 1.0-
i=\
i=\
Eq4.4
where y and y refer to measured and simulated output variables, respectively, and .y is the
1 N
average of the TV measured values (y = —^yt).
It should be noted that in the current version of the SWAT model, the in-stream nutrient
processes have not been linked with the model components that pertain to sediment routing
in the channel network. This is problematic, particularly for total P computations, because
phosphorus is typically bound to sediment and is carried out of most catchments with
sediments. In this study, the results obtained for sediment yield are used for identification of
control processes and critical management actions for control of total P yield.
4.5. Results
Figures 4.2-4.3 depict the cumulative marginal distributions for the sediment and total N
related parameters with highest rank based on the Kolmogorov-Smirnov (dm,n) statistic for
the system behavior. The dm,n test statistic was determined from largest vertical separation
between the cumulative marginal distribution of behavior and non-behavior parameter
distributions. Kolmogorov-Smirnov test statistic (dm,n) has been graphically illustrated on the
first panel in Figure 4.2. Tables 4.1-4.4 summarize the results of the RSA procedure for
sensitivity of input factors listed in Table 2.2. Figures 4.4-4.7 depict the TDSE diagram for
sediment and total N at the outlet of the study watersheds, while Tables 4.5-4.8 provide a
summary of the diagrams.
40
-------
The top three ranks for sediment yield at the outlet of Dreisbach and Smith Fry watersheds
were channel Manning's number CH-N2, peak rate adjustment factor PRF, and average main
channel slope CH-S2. The ranking indicated that fluvial processes within the main channel of
the study watersheds were the control processes for meeting the specified behaviors. Main
channel sediment processes within the SWAT model include channel deposition and channel
erosion (degradation). Whether a channel segment is undergoing erosion or deposition is
determined by comparing its estimated transport capacity with the sediment concentration in
the streamflow. The transport capacity of channel segments is estimated as a function of peak
flow velocity. SWAT uses Manning's equation for estimation of flow velocity in the channel
network. Then, a peak rate adjustment factor is applied to determine the peak velocity in the
segment. Channel Manning's number CH-N2, peak rate adjustment factor PRF, and average
slope of the channel segment CH-S2 interact directly in the mathematical expression of peak
flow velocity of channel segments. Such formulation allows for a negative correlation
between channel Manning's number and peak rate adjustment factor. This negative
correlation can be graphically identified by interpretation of their spider-plots in Figure 4.2.
Sediment yield at the outlet of the study watersheds decreased with increasing CH-N2 value,
while an inverse trend was observed for PRF. In the case of directly correlated model
parameters, without desire for their true values, it is sufficient to adjust only one of the
interacting parameters, and leave the others fixed. Nevertheless, the relative importance of
these parameters leads to the same inferences about the control processes and management
actions. Indirect correlations occur among state variables of the model such as surface runoff
and sediment yield at the outlet. The RSA procedure is essentially a univariate analysis of
system behavior and only suggests the main effects of input factors. Therefore, it is incapable
of revealing any form of correlation structure among directly or indirectly correlated
parameters. A multivariate analysis such as the TSDE procedure was utilized for such
purposes.
Figures 4.4-4.5 show the TSDE trees for sediment yield at the outlet of Dreisbach and Smith
Fry for system behavior specified by the observed past. The summary of main attributes of
the terminal nodes is provided in Tables 4.5-4.6, sorted in order of descending relative
density. High Density Terminal Nodes (HDTN) were selected from the top of the list such
that they contained at least 75% of the input parameter sets (points); i.e. the sum of their
41
-------
corresponding percentage volumes was more than 75 (%). All three HDTNs for the
Dreisbach watershed were defined by the sediment transport process parameters {CH-N2;
PRF; CH-S2}. The factors related to percolation {SOL-AWC} and surface runoff {CN2}
processes defined two HDTNs, but at a lower level in the tree. Likewise the sediment
transport process parameters {CH-N2; PRF; CH-S2} defined both HDTNs identified for the
Smith Fry watershed, while the percolation process parameter {SOL-AWC} defined only one
of them at a lower level. All of these input factors appeared in the RSA classification tables
(Tables 4.1-4.4), except for the available water capacity of the soil layer (SOL-AWC) for the
observed-past behavior defined for Smith Fry. Thus, the multivariate TSDE analysis
confirmed the rankings obtained from the univariate RSA method. Based on the RSA-TSDE
procedure, it can be concluded that sediment transport processes were the control processes
for sediment yield at the outlet of the study watersheds. Therefore, implementation of Best
Management Practices (BMPs) that influence fluvial processes within the channel network
such as grassed waterways and grade stabilization structures would be appropriate choices to
reduce sediment yield at the outlet of the Dreisbach and Smith Fry watersheds. The same
conclusions could possibly be drawn for total P yield because phosphorus is mainly
transported as a sediment bound constituent. A similar analysis can be performed for total P,
once the computation of phosphorus transport through the channel network in SWAT is
linked with sediment transport components of the model.
Average slope steepness of upland fields, parameter SLOPE, had the top RSA rank for total
N load for observed-past behavior in the Dreisbach watershed (see Table 4.5). Variable
SLOPE is a representative topographic attribute of upland areas in the watershed, and is used
by SWAT to estimate the USLE topographic factor used for computation of sheet erosion
from upland areas by Modified Universal Soil Estimation Equation (MUSLE). Other factors
in the MUSLE equation including USLE support practice factor and USLE cover factor also
appear in the list of top RSA ranks. Therefore, sheet erosion from upland areas could be
deemed as the critical process that controls total N yield at the outlet of the Dreisbach
watershed for the specified behavior. Another high ranking parameter was the initial
concentration of organic nitrogen in the soil layer in agricultural areas, denoted by ORGN-
AG. This indicated that the nitrogen loads generated at the agricultural fields were important
to match system behavior. Although ORGN-AG refers to the initial concentration of N in the
42
-------
soil layer, its importance reveals the prominence of anthropogenic agricultural activities such
as fertilizer and manure application, which increase the amount of nitrogen susceptible to
being transported by surface runoff in upland areas. Figure 4.6 depicts the TSDE tree for
total N yield at the outlet of Dreisbach for the observed-past behavior, while Table 4.7
summarizes the main attributes of the terminal nodes and three identifiable HDTNs. The
topographic attribute of upland areas represented by parameter {SLOPE} defined all
HDTNs, suggesting that implementation of BMPs such as parallel terraces that reduce upland
slope steepness would be critical management actions for control of total nitrogen yield at the
outlet of Dreisbach watershed. The other sheet erosion process parameters {USLE-K}, and
also parameters representing agricultural activities that increase concentration of nitrogen in
the soils such as fertilizer application {ORGN-AG} defined two HDTNs in a lower tree level.
Consolidation of the RSA and TSDE results suggest that management actions that focus on
reduction of nitrogen loadings from upland areas are critical to meet the desired target values
in the Dreisbach watershed.
The results of the TSDE procedure for total N yield at the outlet of the Smith Fry watershed
are shown in Figures 4.7 and Table 4.8. Suggestions similar to those made for the Dreisbach
watershed can be inferred for the Smith Fry watershed.
4.6. Discussion
Consolidation of the RSA-TSDE results shows great promise for identification of critical
processes and key management actions through analysis of uncertainty of model simulations.
Particularly, evaluation of effectiveness of non-point source pollution control scenarios such
as implementation of Best Management Practices (BMPs) could be achieved by applying this
method. The common modeling approach is to rely on outputs of a calibrated process-based
model to investigate the impact of management scenarios on pollutant transport within
watersheds. However, this approach can be criticized as a result of its inability to identify a
unique set of process parameters that describe the behavior of the system, often referred to as
problems of identifiability and multiple optima. The success of the combined RSA-TSDE
method is in its ability to identify the critical processes and management actions without
being handicapped by limitations of common parameter estimation procedures. Model
43
-------
realizations obtained from the Monte Carlo routine in the RSA-TSDE procedure can be
further utilized for calculating the margin of safety (MOS) in development of total daily
maximum loads (TMDLs) as suggested in NRC (2001) and Benaman and Shoemaker (2004).
4.7. Conclusions
The conventional utility of process-based watershed models (calibrate —>• validate —>•
predict} for evaluating the impact of management scenarios on fate and transport of
sediments and nutrients was demonstrated in previous sections. In this section, a
computational framework was developed in which investigation of uncertainty provides
complementary quantitative and qualitative information in support of management and
decision making. The analysis focused specifically on two issues. First, the univariate
Regionalized Sensitivity Analysis (RSA) ranking method was applied in conjunction with the
multivariate Tree Structured Density Estimation (TSDE) analysis to identify critical
hydrologic and water quality processes that control the ability to attain specified water
quality conditions. The implications of such analysis provide benefits for watershed
management programs such as TMDL. Under TMDL agenda, Best Management Practices
(BMPs) should be implemented to reduce pollutant loads into the water bodies. Because
BMPs are typically implemented under a restricted budget, identification of critical processes
and management actions becomes highly desirable. Application of the RSA-TSDE procedure
for two small watersheds in Indiana, Dreisbach (6.25 km2) and Smith Fry (7.3 km2), showed
that fluvial channel processes appeared to control sediment yield at the outlet of the
watersheds. Therefore, implementation of within-channel BMPs such as grassed waterways
and grade stabilization structures are likely more effective for reducing sediment loads at the
outlet of the study watersheds. However, total nitrogen yield at the outlets was controlled by
upland nitrogen loading, indicating that implementation of parallel terraces or fertilizer
application strategies would be more effective. Also, EPA has been advised to promote
research on translating effectiveness needs into Use Attainability Analysis (UAA)
implications (NRC, 2001). The adoption of UAA, which is applied for setting new water
quality standards or revising the existing ones, has been restricted due to inadequate technical
guidance from EPA. The results of the TSDE method can provide a numerical probability
44
-------
value that can be utilized to answer pragmatic questions such as how feasible a target value is
for a desired constituent.
The developed methodology can be utilized to inform management about attainability of
desired goals for a given watershed, and key management action(s) for sediment and nutrient
non-point source control. Optimization techniques should be applied to determine the [near]
optimal selection of these management actions based on cost-benefit analysis. For example,
the results of the RS A-TSDE procedure indicated that within-channel BMPs such as grassed
waterways and grade stabilization structures would be most effective for sediment control in
both Dreisbach and Smith Fry watersheds. On the other hand, implementation of upland
nitrogen control plans such as parallel terraces and field borders was deemed more effective
for total nitrogen reduction.
45
-------
Dreisbach Watershed
Smith Fry Watershed
0.2 0.4 0.6 0.8
Normalized CH-N2 Value
Dreisbach Watershed
0.2 0.4 0.6 0.8
Normalized PRF Value
Dreisbach Watershed
0.6 0.8
Normalized CH-S2 Value
Dreisbach Watershed
0.2 0.4 0.6 0.8
Normalized CH-N2 Value
Smith Fry Watershed
0.2 0.4 0.6 0.8
Normalized PRF Value
Smith Fry Watershed
0.2 0.4 0.6 0.8 1
Normalized USLE-K Value
Smith Fry Watershed
0.2 0.4 0.6 0.8
Normalized SOL-AWC Value
0.2 0.4 0.6 0.8
Normalized USLE-P Value
—»- Posterior Behavior Distribution --Q- Posterior Non-behavior Distribution Uniform Distribution
Figure 4.2 Posterior distributions of behavior and non-behavior sets of input factors
along with their prior distribution (uniform distribution) for sediment at the outlet of
Dreisbach and Smith Fry watersheds.
46
-------
Dreisbach Watershed
Smith Fry Watershed
1
c
o
I °-8
f 0.6
5 0.4
1 0.2
0.2 0.4 0.6 0.8
Normalized SLOPE Value
Dreisbach Watershed
0.2 0.4 0.6 0.8
ormalized ORGN.^ Value
AG i
Dreisbach Watershed
0.2 0.4 0.6 0.8
Normalized CN2 Value
Dreisbach Watershed
0.4 0.6 0.8 1
Smith Fry Watershed
0.2 0.4 0.6 0.8 1
Normalized USLE-K Value
Smith Fry Watershed
0.2 0.4 0.6 0.8 1
Normalized USLE-P Value
Smith Fry Watershed
0.2 0.4 0.6 0.8
Normalized USLE-K Value
0.2 0.4 0.6 0.8 1
Normalized USLE-C Value
—*- Posterior Behavior Distribution -o- Posterior Non-behavior Distribution Uniform Distribution
Figure 4.3 Posterior distributions of behavior and non-behavior sets of input factors
along with their prior distribution (uniform distribution) for total N at the outlet of
Dreisbach and Smith Fry watersheds.
47
-------
1
S2
2.249
PRF
I
0.467
PRF
1
S7
0.977
CH-N2
/"sT
[ 0.123
V 5.18
1 J-.
Sn
2.461
SOL-AWC
fSw
/ 0.328
V 1.10
LEGEND
CH Intermediate node
O Low density terminal node
High density terminal node
Figure 4.4 TSDE diagram for sediment at the outlet of Dreisbach watershed:
->nd
rd
1s line-terminal number, 2" line-relative density of the node, 3r line-input factor that splits
the intermediate nodes, or percentage volume of terminal nodes.
S2
1.025
PRF
S3
0.303
PRF
LEGEND
CH Intermediate node
O Low density terminal node
High density terminal node
Figure 4.5 TSDE diagram for sediment at the outlet of Smith Fry watershed, description of
values in each line as in Figure 4.4.
48
-------
1
S5
0.50
USLE-K
XsT
I 1.30
V 23.1
LEGEND
CH Intermediate node
O Low density terminal node
'High density terminal node
Figure 4.6 TSDE diagram for total N at the outlet of Dreisbach watershed, description of
values in each line as in Figure 4.4.
LEGEND
CH Intermediate node
OLow density terminal node
'High density terminal node
Figure 4.7 TSDE diagram for total N at the outlet of Smith Fry based on observed-past
behavior definition, description of values in each line as in Figure 4.4.
49
-------
Table 4.1 Regionalized Sensitivity Analysis ranking of input factors for sediment based on
observed-past behavior definition at the outlet of Dreisbach watershed
Symbol
Short Description
CH-N2 Channel Manning's number
PRF Peak rate adjustment factor
CH-S2 Average main channel slope
SOL-A WC Available water capacity of soil layer
CN2 SCS curve number
ESCO Soil evaporation compensation factor
SPEXP Exponent coefficient for sediment routing
ALPHA-BF Baseflow alpha factor
Summary
Statistics:
Number of Behavior Simulations (m) :
Total number of Simulation (m+n) :
Success rate :
Rank
&m,n
0.4822
0.3090
0.2189
0.1457
0.1346
0.1056
0.0730
0.0532
1010
5000
20.2%
Table 4.2 Regionalized Sensitivity Analysis ranking of input factors for sediment based on
observed-past behavior definition at the outlet of Smith Fry watershed
Symbol
Short Description
Rank
CH-N2 Channel Manning's number
PRF Peak rate adjustment factor
CH-S2 Average main channel slope
USLE-K USLE soil erodibility factor
USLE-P USLE support practice factor
0.5929
0.3616
0.2010
0.1172
0.0970
Summary Statistics:
Number of Behavior Simulations (m) : 1360
Total number of Simulation (m+n) : 5000
Success rate : 27.2%
50
-------
Table 4.3 Regionalized Sensitivity Analysis ranking of input factors for total N based on
observed-past behavior definition at the outlet of Dreisbach watershed
Symbol
SLOPE
ORGNAG
CN2
USLE-K
SOL-AWC
USLE-P
ORGNPAST
SLSUBBSN
USLE-C
Short Description
Average slope steepness
Initial organic N in soil, agricultural land use
SCS curve number
USLE soil erodibility factor
Available water capacity of soil layer
USLE support practice factor
Initial organic N in soil, pasture land use
Average slope length
Maximum value of USLE cover factor
Summary Statistics:
Number of Behavior Simulations (ni) :
Total number of Simulation (m+n) :
Success rate :
Rank
Ufn^n
0.2353
0.1755
0.1539
0.1345
0.1266
0.1096
0.0877
0.0593
0.0550
1448
5000
29.0%
Table 4.4 Regionalized Sensitivity Analysis ranking of input factors for total N based on
observed-past behavior definition at the outlet of Smith Fry watershed
Symbol
Short Description
Rank
ORGN AG Initial organic N in soil, agricultural land use
USLE-K USLE soil erodibility factor
USLE-P USLE support practice factor
USLE-C Maximum value of USLE cover factor
CH-K2 Channel effective hydraulic conductivity
0.4436
0.3848
0.2514
0.0908
0.0518
Summary Statistics:
Number of Behavior Simulations (m) : 2039
Total number of Simulation (m+n) : 5000
Success rate : 40.8%
51
-------
Table 4.5 Summary of TSDE terminal nodes for sediment related parameters based on
observed-past behavior definition at the outlet of Dreisbach watershed
Terminal Relative
Node Density 1
$5
814
S4
89
Sw
86
High
5.719
2.903
2.451
1.491
0.723
0.414
0.328
5.18
CH-N2
CH-N2
CH-N2
CH-N2
CH-N2
CH-N2
CH-N2
CH-N2
Density Terminal Nodes
Tree level (root node =
234
PRF
PRF
PRF
PRF
PRF
PRF
PRF
PRF
(HDTNs): {
CH-N2 CH-S2
CH-N2 CH-S2
CH-N2 CH-S2
CH-N2
CH-N2 CH-S2
S^S.-Su-S^}
level 1)
5 6
SOL-AWC CN2
SOL-AWC CN2
SOL-AWC
Table 4.6 Summary of TSDE terminal nodes for sediment related parameters based on
observed-past behavior definition at the outlet of Smith Fry watershed
Terminal Relative
Node Density i
S5 1.90 CH-N2
Sn 1.40 CH-N2
Sw 0.942 CH-N2
S4 0.64 CH-N2
S9 0.301 CH-N2
S6 0.125 CH-N2
iigh Density Terminal Nodes
2
PRF
PRF
PRF
PRF
PRF
PRF
(HDTNs):
Tree Level (root node = level 1)
3 4
CH-S2 SOL-AWC
CH-S2 SOL-AWC
CH-S2
{S5.Sn}
52
-------
Table 4.7 Summary of TSDE terminal nodes for total N related parameters based on
observed-past behavior definition at the outlet of Dreisbach watershed
Terminal Relative
Node Density 1
$2
89
High
1.90
1.40
1.30
0.90
0.3
SLOPE
SLOPE
SLOPE
SLOPE
SLOPE
Density Terminal Nodes
Tree level (root node = level 1)
234
ORGNAG USLE-K
ORGNAG
ORGNAG USLE-K
ORGNAG USLE-K
(HDTNs): {S^S2}
SLOPE
SLOPE
Table 4.8 Summary of TSDE terminal nodes for total N related parameters based on
observed-past behavior definition at the outlet of Smith Fry watershed
Terminal Relative
Node Density 1
Tree level (root node = level 1)
1.836 ORGNAG USLE-K
0.410 ORGNAG USLE-K
0.231 ORGNAG
High Density Terminal Nodes (HDTNs): { S5}
53
-------
SECTION 5. COST-EFFECTIVE ALLOCATION OF WATERSHED MANAGEMENT PRACTICES
USING A GENETIC ALGORITHM
5.1. Abstract
Implementation of conservation programs are perceived as being crucial for restoring and
protecting waters and watersheds from nonpoint source pollution. Success of these programs
depends to a great extent on planning tools that can assist the watershed management
process. Herein, a novel optimization methodology is presented for deriving watershed-scale
sediment and nutrient control plans that incorporate multiple, and often conflicting,
objectives. The method combines the use of a watershed model (SWAT), representation of
best management practices, an economic component, and a genetic algorithm-based spatial
search procedure. For two small watersheds in Indiana located in the Midwestern portion of
the United States, selection and placement of best management practices by optimization was
found to be nearly three times more cost-effective than targeting strategies for the same level
of protection specified in terms of maximum monthly sediment, phosphorus, and nitrogen
loads. Conversely, for the same cost, the optimization plan reduced the maximum monthly
loads by a factor of two when compared to the targeting plan. The optimization methodology
developed in this study can facilitate attaining water quality goals at significantly lower costs
than commonly used cost-share and targeting strategies.
5.2. Introduction
Best management practices (BMPs) are widely accepted as effective control measures for
agricultural nonpoint sources of sediments and nutrients. The 2002 Farm Bill provided up to
$13 billion for conservation programs aimed at protecting water quality from agricultural
nonpoint source (NFS) pollution (USDA, 2003). In addition, under the Clean Water Act
Section 319 Nonpoint Source National Monitoring Program and wetland protection
programs, the EPA supports programs to reduce the negative impacts of runoff from
agricultural, urban, and industrialized areas. Similarly, the Natural Resources Conservation
54
-------
Service (NRCS) provides hundreds of millions of dollars in federal funds to support
agricultural best management practices (BMPs) in an effort to reduce the movement of
pollutants into our waterways. Success of such programs, however, is contingent upon
availability of efficient watershed-scale planning tools.
Implementation of BMPs is challenged by complexities in incorporation of conflicting
environmental, economic, and institutional criteria. Environmental assessments in watersheds
hinge on resolving social benefits such as achieving the goal of swimable and fishable water
bodies under the EPA's Total Maximum Daily Load (TMDL) agenda. While BMPs facilitate
achievement of such targets, their establishment bears additional cost for watershed
management and/or agricultural producers. Since management practices are usually
implemented under a limited budget, costs associated with unnecessary/redundant
management actions may jeopardize attainability of designated water quality goals.
Identifying optimal combinations of watershed management practices requires systematic
approaches that allow decision makers to quickly assess trade-offs among environmental and
economic criteria.
Cost-sharing with landowners is promoted by government agencies for BMP implementation
in agricultural fields (EPA, 2003). BMPs are implemented through site investigation,
monitoring, and field-scale modeling. While cost-share programs may improve water quality
standards at a field scale, their impact at a watershed scale typically remains unknown for
two main reasons. First, impact of BMPs may duplicate or overlap each other, thereby
reducing the potential benefit for the watershed. Interactions between BMPs may also
significantly affect their individual performances at a watershed scale. Second, water quality
impacts of BMPs are site-specific, greatly influenced by landscape characteristics such as
land use, soils, and management actions. Thus, the same BMP is likely to have varying
efficacies at different locations within a watershed. The direction of many regulatory
programs such as the EPA's Clean Water Act (CWA) has now shifted from a source-by-
source, pollutant-by pollutant approach to a more holistic watershed-scale strategy (EPA,
2003).
55
-------
An example of watershed scale approaches is targeting where pollution control measures are
allocated to critical source areas. Critical sources are portions of the watershed that are
believed to contribute intensively to nonpoint source pollution. Previous studies in the same
watershed, watershed data (i.e. land use, soils, and management characteristics), physical
characteristics such as topography and proximity to a stream, and instream monitoring data
are typically analyzed to identify critical sources. This approach has been endorsed in the
recent draft of EPA's "Handbook for developing watershed plans to restore and protect our
waters" (EPA, 2005a). However, identification of an optimal pollution control strategy
through targeting becomes infeasible in large complex watersheds, because the number of
possible scenarios increases exponentially with the number of fields. The search for an
optimal watershed pollution control plan needs to be conducted in a more efficient manner.
Development of increasingly powerful computers has facilitated application of mathematical
programming heuristics for solving complex and computationally cumbersome problems. A
few previous studies have used evolutionary algorithms (EA) to optimize placement of
watershed management actions. Srivastava et al. (2002) linked a genetic algorithm (GA) with
a continuous NPS model (AnnAGNPS) to optimize selection of crop rotation practices in a
7.25 km2 USDA experimental watershed in Pennsylvania. The optimized crop management
scheme obtained after 150 GA generations decreased pollution load by nearly 56% and
increased net annual return by nearly 110% from a random selection of crop rotation
scenarios. The objective function of the optimization procedure was formulated such that
either pollution reduction or net return, but not both at the same time, was maximized. Veith
et al. (2004) developed a GA-based optimization model for cost-effective allocation of land
use and tillage practices to upland fields that minimized cost (economic criteria) while
sediment load at the outlet of watershed was constrained to a predefined value. The model
was tested in a 10.14 km2 case study watershed in Virginia. Results of the study showed that
the optimization plan achieved the same sediment reduction as a targeting plan at a lower
cost. The major setback with optimizing only one of the decision objectives (economic,
environmental, and/or institutional) is that the procedure can not demonstrate the tradeoff
between multiple objectives.
56
-------
Incorporation of multiple objectives in the search for alternative agricultural landscapes was
studied by Bekele and Nicklow (2005) and Muleta and Nicklow (2005). A mutiobjective
evolutionary algorithm (SPEA; Strength Pareto Evolutionary Algorithm, Zitzler and Thiele.,
1999) was linked with a watershed model for selection of crop rotation and tillage operations
in a 133 km2 watershed in southern Illinois. Multiple objectives included minimizing
pollutant loads, and maximizing net profit. Though both studies showed the effectiveness of
an integrated approach in providing tradeoff between multiple objectives, explicit
incorporation of water quality and budget constraints was not examined. These two studies
do not offer any clear guidance for selection of operating parameters and termination criteria
as they rely on the previously developed SPEA. Moreover, neither optimization results were
compared with the more commonly used targeting strategies, nor was convergence of
pollutant loads/net profit tested. Hence, appraisal of efficiency of the optimization method in
deriving solutions better than targeting strategies was curtailed.
Research to date indicates the promise of heuristic optimization for cost-effective allocation
of watershed management practices. Unlike gradient-based approaches, heuristic techniques
do not require linearity, continuity, or differentiability either for objective/constraint
functions or for input parameters. Thus, they are well-suited for cost-effective allocation of
watershed management plans. However, several questions still defy answers. A decision
making tool that can clearly accommodate economic, environmental and institutional criteria
is still lacking. The means for imposing target values for pollutant loads, and total watershed
cost of implementation of management plans needs to be explored. As discussed by
Srivastava et al. (2002), a more versatile formulation is needed to maximize pollution
reduction and minimize cost at the same time. Moreover, previous studies have only focused
on alternative agricultural landscapes. Research on economic evaluation of structural BMPs
that are installed both in upland areas and within the channel network is lacking. In such
cases, the interaction between the BMPs could be critical to reach water quality goals.
Finally, selection of GA's operating parameters, and termination criteria need further study.
The main goal of this section is to develop an optimization framework that enhances decision
makers' capacity to evaluate a range of agricultural and environmental management
alternatives. The tool will be designed to identify near optimal watershed plans that reduce
57
-------
pollutant loads at a watershed outlet to below regulatory or target values with minimum cost.
We hypothesize that reductions of pollutants at watershed outlets can be attained at
significantly lower cost by optimized implementation of conservation practices than by
current cost-share and targeting approaches. This overall goal is achieved by the following
specific steps:
1. Development of a novel genetic algorithm-based spatial search model. This step will
focus on formulating versatile objective and constraint functions for the optimization
model that can handle multi-criteria and landscape characteristics.
2. Integration of an NFS model (SWAT; Soil and Water Assessment Tool), a new BMP
representation method, and a cost-benefit economic relationship with the GA-based
spatial optimization model to identify optimal spatial allocation of best management
practices;
3. Provision of guidance for selection of GA operating parameters and termination criteria;
4. Conducting case studies to determine how the optimization plan compares to plans from
cost-share and targeting strategies.
5.3. Theoretical Considerations
5.3.1. Genetic Algorithm (GA)
Genetic algorithms belong to the evolutionary class of artificial intelligence (AI) techniques.
GAs are based on natural selection of chromosomes from a population for mating,
reproduction of offspring by crossover, and mutation to ascertain diversity- ideas borrowed
from biology. Each chromosome string in the population corresponds to a solution for the
problem at hand, with each variable being represented by a gene (a specific position in the
string). The values of the genes, known as allele, can be binary, real-valued, or character-
valued.
The GA begins with an initial population containing randomly generated chromosomes or
strings, each representing a possible solution to the problem. The next population is
generated by mating the fittest solutions (i.e. chromosomes) in the previous population.
58
-------
Based on the principle of survival of the fittest, the higher the solution's fitness, the more
likely it will be chosen for reproduction. The search for the optimal solution continues until a
predefined termination criterion is reached.
The GA is essentially a search algorithm for solving difficult nonlinear optimization
problems and is not intended to examine all possible solutions. Regardless of how long the
algorithm searches, its convergence to an optimum can not be guaranteed. However, it has
been shown that the GA converges to near optimal solutions for a variety of problems
(Winston and Venkataramanan, 2003). Efficiency of the algorithm depends on the
optimization's operating parameters and the convergence criterion. The higher the number of
individual evaluations for converging to the optimum, the less efficient is the procedure. The
values of the operating parameters are problem dependent, and can be determined by
performing a sensitivity analysis.
The GA is a well-suited optimization technique for spatial allocation of BMPs, because
unlike gradient-based methods it does not require linearity, continuity, or differentiability of
either the objective function or the constraint function.
5.4. Methodology
The optimization model developed in this study is comprised of the SWAT model for
simulating pollutant loads, a BMP representation tool, an economic component, and a GA-
based spatial optimization technique. A MATLAB (The Mathworks, Natick, MA) computer
program was developed to provide the linkage among various components of the model as
shown in Figure 5.1. The model was tested for optimization of the location of field borders,
parallel terraces, grassed waterways, and grade stabilization structures in the Dreisbach and
Smith Fry watersheds.
SWAT simulations were performed for a 10-year period from January 1st, 2000 through
December 31st, 2009. In the analysis, 1991-2000 precipitation data, 2000 USDA-National
Agriculture Statistics Service (NASS) land use, and 2002 Soil Survey Geographical Database
(SSURGO) were utilized to establish a base-case SWAT run. Parameter values in the base-
case run were selected from a manual calibration (Arabi et al., 2005). Portions of the
59
-------
Start here
First generation:
Choose randomly
Nest generations:
Crossover/
Replacement/
Mutation
ergence of fit
and constraints
Water duality tareets
On-site sediment and
nutrient reduction
benefits
Cost-benefit ratio
Establishment,
maintenance, and
opportunity cost
Figure 5.1 Schematic of the optimization procedure.
watershed classified as urban and forested areas were not considered for implementation of
BMPs.
5.4.1. NFS Model
The optimization tool developed here does not necessarily require using the SWAT model as
the NFS component. The SWAT model was chosen because the case studies will deal with
sediments, and nutrients. Soundness of mathematical representation of sediment, nutrient,
and pesticide processes in SWAT has been validated in previous research (Santhi et al., 2001;
Kirsch et al., 2002; Arabi, 2005; Vazquez-Amabile, 2005). Moreover, SWAT has been
linked with optimization routines for automated calibration of the model (van Griensven and
Bauwens, 2003) and for evaluation of efficiency of nonpoint source pollution regulatory
programs (Whittaker et al., 2003). The optimization framework herein can be easily
integrated with other hydrologic/water quality models to develop management plans for other
types of pollutants, and for other watersheds as well.
60
-------
5.4.2. BMP Representation
In this study, a method presented by Arabi et al. (2004) and Bracmort et al. (2006) was
utilized to evaluate the water quality impacts of grassed waterways, grade stabilization
structures, field borders and parallel terraces. The method was developed based on published
literature pertaining to BMP simulation in hydrological models and considering the
hydrologic and water quality processes simulated in SWAT. Based on the function of the
BMPs and hydrologic and water quality processes that are modified by their implementation,
corresponding SWAT parameters were selected and altered. Table 2.3 summarizes the
SWAT parameters and their corresponding values for representation of the BMPs. Arabi et
al. (2004) provides a detailed description of the method used for representation of field
borders, parallel teraces, grassed waterways, and grade stabilization structures.
5.4.3. Economic Component
An economic component was developed for the optimization model that is comprised of a
cost function in addition to an economic return (benefit) function. Both cost and benefit are a
function of watershed characteristics and time. The benefit function reflects the impact of
BMPs on sediment and nutrient reductions.
Cost Function
The total cost of implementation of BMPs was evaluated by establishment, maintenance, and
opportunity costs. Establishment costs included the cost of BMP installation, and technical
and field assistance. Maintenance cost is usually evaluated as a percentage of establishment
cost. The opportunity cost is a dollar value that would be produced over the BMP design life
as a result of investing the establishment and maintenance costs by purchasing saving bonds.
For each individual BMP, the total cost (cf) was evaluated by the following equation:
s-} EqS.l
T=\
c0=c0(x,0 Eq5.2
61
-------
where c0 is the establishment cost that is a function of state of the watershed (i.e. x, landscape
and physical characteristics) at time t; s is the interest rate; td is BMP design life; and rm is
the ratio of maintenance cost to establishment cost. Interest rate (s) of 6.5% was used in
computations. The design life of a BMP is defined by Natural Resources Conservation
Service (NRCS) as "the intended period of time that the practice will function successfully
with only routine maintenance determined during design phase" (USDA-NRCS, 2005).
Design life, establishment cost, and maintenance rate for BMPs in this study are summarized
in Table 5.1. The table also provides information with regard to the number and area (i.e.
quantity of each type of BMP) allocated within the study watersheds through targeting (see
Figure 2.1). For each management plan, establishment cost (c0) of each type of BMP was
obtained by multiplying its unit establishment cost (c) by its quantity (abmp). The total
watershed cost of BMPs was computed by summing up individual costs from Eq 5.1.
Table 5.1 Unit cost of establishment (c), maintenance rate (rm), and design life (td) of
BMPs in the study along with number (nbmp) and quantity (abmp) of each type in the study
watersheds as shown in Figure 2.1 (i.e. targeting strategy).
BMP
Field Border
Parallel Terrace
Grassed Waterway
Grade Stabilization
a
($)
4.6 /m
26 /m
6200 /ha
10,000 /structure
rm -
1
3
3
2
, Dreisbach
(yr)
10
10
10
15
nbmp
7
4
5
10
Smith Fry
abmp ,
, .\. nbmp
(unit) ^
2600 (m)
2130 (m)
3. 50 (ha)
1
2
1
2
abmp
(unit)
1800(m)
480 (m)
0.95 (ha)
- Adapted from Indiana Environmental Quality Incentives Program (2004).
Benefit Function
The economic return of implementation of BMPs was determined by assigning monetary
values to onsite and offsite benefits of sediment and nutrient reductions. Reduction of
sediments and nutrients as a result of implementation of a management practice is a function
of landscape characteristics. Likewise, benefits gained by BMP implementation vary with
landscape properties of the site where the BMP is installed. While offsite benefit of BMPs
were determined based on cost of offsite damage that would be caused by sediment and
nutrient loads, onsite benefits were estimated based on the quantity of reduced nutrient loads
from fields under the influence of BMPs.
62
-------
Ribaudo et al. (1989) equated the benefit of reducing soil loss to offsite damage in $ per ton
of eroded soil based on freshwater recreation, marine recreation, water storage, navigation
flooding, roadside ditches, irrigation ditches, freshwater commercial fishing, marine
commercial fishing, municipal water treatment, and municipal and industrial use. A $1.15/t
damage cost estimate was suggested for the Corn-Belt region in the United States.
Additionally, Cangelosi et al. (2001) evaluated the benefit of reducing soil erosion for
dredging in the Maumee River basin to be $0.87 for a ton of eroded soil. In a study by Fang
and Easter (2003), a social cost of $2.65/kg of removed phosphorous was used to quantify
the cost-effectiveness of a water quality trade in southeastern Minnesota. Onsite water
quality benefit of BMP implementation was evaluated by the monetary value of nutrient
reduction from upland areas. Buckner (2001) estimated the benefit of reducing phosphorus
from agricultural lands by the bulk rate cost of triple super phosphate (00-46-00) that, in
2003, was $0.26/kg (Heartland Co-op, IN, personal communication). Table 5.2 summarizes
the monetary benefits (bm) of unit reduction of sediment, phosphorus, and nitrogen deliveries
expressed in 2000 dollars used for the Dreisbach and Smith Fry watersheds. The values used
in this study are site specific and may vary for other regions.
After on-site and off-site benefits of BMPs were quantified by monetary values from related
literature, the total benefit (hi) was calculated for individual BMPs over their design life as:
^LV •".'
-------
Table 5.2 Monetary benefits of unit reduction of
sediment and nutrient loads in 2000 dollars
Variable
Sediment
Phosphorus
Nitrogen
Benefit (bni)
($)
2.5 /t
2.86 /kg
3.5 /kg
It is worthwhile noting that bmt serves as a weighting factor that incorporates the relative
importance of reduction of pollutant constituent / in development of watershed management
plans. A higher value indicates that larger benefits can be gained from unit reduction of the
corresponding constituent. In the absence of regional and site-specific data for the sort of
analysis that was used in this study, the term bm can be assigned by judgment of the analyst.
If all constituents bear the same importance, bm should be 1 for all constituents of concern.
5.4.4. Optimization Component
A genetic algorithm (GA) was employed to optimize the spatial allocation of BMPs. In this
GA component, each optimization string (vectora') corresponds to a specific watershed
management plan. Figure 5.2 is a demonstration of an optimization string for placement of
BMPs. The length of each string (m) corresponds to the total number of genes, i.e., individual
management actions (»},/) that are considered in the optimization. For example, in a
watershed with 50 fields (i.e., nhm=5Q) considered for implementation of field borders
and/or parallel terraces, and 20 reach segments (i.e., nch=20) considered for implementation
of grassed waterways and/or grade stabilization structures, the total number of genes on each
management string is equal to m= [2x50+2x20=] 140. The alleles are binary values, with
"1" or "zero" indicating that the corresponding BMP "be" or "not be" implemented. The /'th
optimization generation (f) with nstr solution vectors is expressed as:
P' =
Eq5.5
_nstr
a
nstr,l
64
-------
A string (chromosome)
gss gww pt fb
nch — nch —— nhm —— nhru
Figure 5.2 Schematic of an optimization string (chromosome) representing grade
stabilization structures (gss), grassed waterways (gww), parallel terraces (pt), and field
borders (fb). nch and nhru refer to the number of channel segments and fields (HRUs) in the
watershed, respectively.
where
V/e{l,2,...,7wfr},V/e{l,2,...,y»}:aj./ =0/1 Eq 5.6
where a] denotes the/h management vector in P\ and al} ^reflects the 7th gene (an individual
management action) ma'}. . For a total of ngen optimization generations, the entire solution
space (A) is obtained from the union of solutions in optimization generations:
Mathematical Formulation
A general multi objective optimization problem can be formally expressed as (Zitzler and
Thiele, 1999):
Minimize/Maximize z(a) = f(a) = (_/[ (a), /2 (or),...,/n (or)) Eq 5.8
subject to
a = (a1,a2,...,am)en Eq 5.9
z = (z1,z2,...,zn)eZ EqS.10
65
-------
where z denotes a vector function / that maps a set of m decision variables a to a set of n
objectives; Q and Z are, respectively, the decision parameter space and objective space (refer
to Ringuest, 1992 for more information).
Several common methods have been previously used for handling multiobjective
optimization problems (Deb, 2001). The most commonly used approach is the weighted sum
approach. In this method, a set of objectives are aggregated into a single objective by
multiplying each objective by a user-defined weight (w). These weights reflect the
importance of each objective in the context of the optimization problem:
Z^tftW Eq5.ll
;=1
The BMP spatial allocation problem was formulated with an aggregate single-objective
based on multiple criteria. The multiple criteria included minimizing implementation cost of
BMPs and maximizing pollutant load reductions. Monetary benefits of unit reduction of
sediment, phosphorus, and nitrogen loads (Table 5.2) were used to determine the weighted
sum of onsite and offsite benefits gained by reducing pollutant loads (bf). Then, the objective
function was evaluated with a benefit (bf) to cost (cf) ratio at the watershed level. Pollutant
loads were constrained to ascertain that sediment and nutrient yields would not exceed
allowable values and/or that total available budget is not exceeded by entailed costs.
Allowable pollutant loads are typically specified by regulatory agencies. Since the available
budget for implementation of BMPs is usually limited, an additional constraint was imposed
to search for solutions that can be implemented with a cost less than the specified budget.
The mathematical representation of the objective function used in this study was:
bt
Maximizez = = , T~l ' ^ Eq5.12
r,l) (
max
subject to water quality constraints:
V T=\
66
-------
F! (x, t, td, a) = ymaxt (x, t, td, y, a)-ytt<0 Eq 5.13
and budget constraints:
T2(x,t,td,a) = ct(x,t,td,a)-wcost<0 Eq5.14
where ymaxt is maximum delivery of pollutant constituent /' after implementation of BMP
combination (a) estimated with SWAT simulations over period td; and ytt is the allowable
load of constituent /'. Variable ct is the total cost of implementation of a estimated by
applying Eq 5.1, and wcost is the total available budget for implementation of management
plans. The denominator in 5.12 was designed such that it will never be zero.
The optimization constraints in 5.13 and 5.14 are typically defined by regulatory and
implementation agencies. For example, allowable sediment and nutrient loads (yt,) may be
obtained from a Total Maximum Daily Load (TMDL) for a given watershed. While ytt and
ymaxt can be expressed on a daily, monthly, or annual basis, as loads or concentrations, their
units should be consistent. The cost constraint represents available budget for implementation
of watershed management scenarios and may be specified by implementation agencies.
The fitness score (fs) for each string was evaluated by the objective function (z) associated
with the string from Eq 5.12. Infeasible solutions, i.e., solutions that do not satisfy the
constraint functions FI and F2 in 5.13 and 5.14, were penalized by applying a penalizing
factor k as:
fs = kxz;
10~5 ^vr^O Eq5.15
k =
II
Operating Parameters
Selection of the size (nstr) of the optimization population (P) and their initial values (cr)y) is
complex and is likely to vary for different problems (Deb, 2001). A large population provides
the GA with an adequate sampling of the search space. However, there are some trade-offs
67
-------
between the population size and the number of generations needed for convergence. A small
population size may cause the GA to become trapped in a local optimum, whereas a large
size may be computationally inefficient and take too long to converge. Winston and
Venkataramanan (2003) suggest that a population of 50 or 100 is typical for GAs. In this
study, population size (nstr) was determined through a sensitivity analysis.
The initial population (Po) was generated by assigning random binary values (i.e. 0 or 1) for
management actions in each watershed plan as shown in Figure 5.2. Subsequent generations
were produced through the following steps:
1 . Replacement: The strings of the previous population characterized by the best fitness
scores were propagated to the next generation unchanged. The rest of the population was
replaced by new chromosomes generated from mating. The replacement rate (rr) was
determined through a sensitivity analysis.
2. Selection: The parents from the previous generation were selected to reproduce the
offspring through the principle of survival of the fittest. A probability score (F) was
assigned to each string based on its fitness score (fs in Eq 5. 15) as:
Eq5.16
where nstr is the population size. Then, the strings were sorted in descending order of the
assigned probabilities, and their cumulative probabilities were determined. Two random
numbers were generated for selection of two parent chromosomes for mating.
3. Mating: The two selected parents were used for production of two offspring. A crossover
point was selected randomly so that the offspring consists of a portion of each parent. In
this study, a random crossover point was chosen for each type of BMP. One point was
chosen to crossover the genes corresponding to parallel terraces, one point for field
borders, one for grassed waterways, and another one for grade stabilization structures.
68
-------
4. Mutation: The genetic makeup of the offspring was altered randomly to avoid local
optima. Typically 1% to 5% of the genes mutate per iteration. A 2% mutation rate (mr)
was used here. The alleles for the randomly chosen genes changed from a "1" to "0", or
vice versa.
5. Repeat: Previous steps were repeated until the algorithm converged.
The optimization procedure was terminated once the convergence criteria, defined below,
were met: (1) the fitness score of the best solution of ten consecutive generations did not
vary; (2) the median fitness score of solutions in 10 consecutive generations did not vary
more than 5%; and (3) values of constraint functions did not vary for ten consecutive
generations. A minimum of 2000 individual model evaluations were considered.
Efficiency of the optimization algorithm for various combinations of operation parameters
was evaluated through a sensitivity analysis. Table 5.3 summarizes different combinations of
operating parameters (labeled as setup 1 to setup 6) that were examined in the analysis. A
total of 2000 individual model evaluations were performed for each combination of operating
parameters. As a result, all setups had almost identical computational runtime. The number of
generations (ngen) in Table 5.3 was determined as:
ngen =
nmax
nstr • rr
Eq5.17
where nmax, nstr, and rr refer to total number of model evaluations, population size, and
replacement rate, respectively. The population size, and replacement rate with the highest
efficiency was utilized to identify the optimal spatial allocation of BMPs in the study
watersheds. The setup with the highest objective function value over all model simulations
was appraised to be the most efficient.
Table 5.3 Various combinations of GA operating parameters for sensitivity analysis
Parameter
Jopulation Size
leplacement Rate (%)
lumber of Generations
Symbol
nstr
rr
ngen
I
100
70
28
2
100
90
22
Setup
3
50
70
58
4
50
90
45
5
20
70
142
6
20
90
110
69
-------
5.4.5. Optimization Cases
BMPs depicted in Figure 2.1 were implemented in the Dreisbach and Smith Fry watersheds
through targeting strategies. This provided the opportunity to compare the cost-effectiveness
of BMP combination obtained from the optimization procedure with the combination
allocated through targeting. Four cases (summarized in Table 5.4) were examined in this
study. Case 1 represented the base-case simulation with no BMPs in place. Cost-
effectiveness of BMPs allocated by targeting was determined in case 2. The results of the
optimization procedure were compared with targeting in two ways. In case 3, the
optimization procedure was used to allocate BMPs such that maximum sediment,
phosphorus, and nitrogen loads over the simulation period (2000-2009) (ymaxt in Eq 5.13)
did not exceed the ones corresponding to the targeting plan. The purpose of comparing cases
2 and 3 was to compare the total watershed cost (ct in Eq 5.14) of the two cases while
providing the same level of water quality protection. For simulations in case 4, the
optimization constraints consisted of only the cost constraint. This case aimed at
investigating environmental benefits of BMPs from optimization and targeting plans limited
to the same budget for implementation. Therefore, the optimization procedure was designed
to terminate once the total watershed cost of the near optimal solution reached the cost of the
targeting plan.
Table 5.4 Description of cases compared in this study
Case Strategy
Description
1 Base-case No BMP implemented in the watershed.
2 Targeting BMPs spatially allocated as in Figure 2.1.
3 Optimization 1 BMPs allocated such that maximum monthly sediment and nutrient
loads did not exceed maximum monthly values from targeting
strategy. No cost constraint was imposed.
4 Optimization 2 BMPs allocated such that total watershed cost did not exceed cost of
targeting strategy.
70
-------
5.5. Results
Results of sensitivity analysis for the GA operating parameters indicated that the algorithm is
more efficient with smaller population size and replacement rate. Figure 5.3 shows sensitivity
of the optimization procedure to combinations of operating parameters listed in Table 5.3.
Model evaluations were performed for the Dreisbach watershed over the 2000-2009 period.
The trend for each setup was determined based on a total number of 2000 individual model
evaluations. The runtime corresponding to 2000 model evaluations for the study watersheds
was nearly 15 hours on a 2.0 GHz PC. The corresponding number of generations for each
setup was determined from Eq 5.17. A mutation rate of 2% was used for all setups. The
results indicate that with the same number of model evaluations and computational cost, the
setup with the highest number of generations (population size of 20 with 70% replacement
rate, i.e. setup 5 in Table 5.3) outperformed the other combinations. Similar results were
obtained for the Smith Fry Watershed. Therefore, these values were used for optimizing the
location of grassed waterways, grade stabilization structures, parallel terraces, and field
borders in the study watersheds.
BMPs allocated with the optimization procedure were more cost-effective than randomly
selected BMPs and BMPs allocated through targeting strategies. Table 5.5 summarizes
results for base-case and targeting cases, case 1 and case 2. Estimated average and maximum
0.4
0.3
(D
^5"
O
0.2
- Setup 1
- Setup 2
- Setup 3
• • Setup 4
- Setup 5
• - Setupd
500 1000 1500
Number of Individual Model Evaluations
2000
Figure 5.3 Analysis of sensitivity of GA operating parameters. Table 5.3 shows values of
operating parameters corresponding to each setup.
71
-------
Table 5.5 Results for base-case and targeting cases. Mean and maximum value of monthly
sediment, phosphorus, and nitrogen loads were estimated from SWAT simulations.
Watershed Variable
Mean
Sediment Max
Reduction
Mean
<-i
g Phosphorus Max
"« Reduction
'(3 -. ,
i? Mean
Nitrogen Max
Reduction
objective function
watershed cost
Mean
Sediment Max
Reduction
Mean
•^_
£ Phosphorus Max
£ Reduction
S Mean
Nitrogen Max
Reduction
objective function
watershed cost
- Used as water quality constraint for case
c , , TT ., Base-case
Symbol Units ,„ ...
J (Case 1)
ys
ymaxs
rs
yP
ymaxp
rp
yn
ymaxn
rn
z
ct
ys
ymaxs
rs
yP
ymaxp
rp
yn
ymaxn
rn
z
ct
t/ha/m
t/ha/m
%
kg/ha/m
kg/ha/m
%
kg/ha/m
kg/ha/m
%
$/$
$
t/ha/m
t/ha/m
%
kg/ha/m
kg/ha/m
%
kg/ha/m
kg/ha/m
%
$/$
$
3 summarized in Table 5
- Used as cost constraint for case 4 summarized in
Table 5.6.
0.074
0.55
-
0.044
0.20
-
0.395
2.69
-
-
-
0.116
1.62
-
0.32
2.129
-
4.545
34.52
-
-
6.
Targeting
(Case 2)
0.022
0.17 a
70
0.033
0.15-
25
0.298
2.1 a
25
0.12
414,690 -
0.105
1.5
10
0.306
2.06
5
4.368
33.31
4
1.4
60,610h
monthly sediment and nutrient yields simulated by SWAT, and total watershed cost (ct) are
provided for the study watersheds. Variable z in the table refers to the benefit-cost ratio as
defined in Eq 5.12. Percent reduction of constituent /' (r,) as a result of implementation of
BMPs was determined as:
r. =
Eq5.18
where /' is the constituent of interest (i.e., sediment, phosphorus , or nitrogen load), yn is the
average monthly load of constituent /' for the base-case (case 1), and>^2 is the average
72
-------
monthly load of constituent /' for the targeting case (case 2). In the Dreisbach watershed,
implementation of the targeting plan (depicted in Figure 2.1) would cost $414,690 over the
2000-2009 period, with a benefit-cost ratio (z in Eq 5.12) of 0.12. As a result, average
monthly sediment, phosphorus, and nitrogen yields would be reduced from the base-case by
nearly 70%, 25%, and 25%, respectively. Likewise, estimated maximum monthly sediment,
phosphorus, and nitrogen yields at the outlet of the Dreisbach watershed would respectively
reduce to 0.17 t/ha, 0.15 kg/ha, and 2.1 kg/ha. These values were used as the allowable
sediment and nutrient loads (ytt in Eq 5.13) for the first optimization case (case 3).
A summary of results for optimization cases in the study watersheds is provided in Table 5.6.
Highlighted values reflect optimization constraints in each case that were obtained from
targeting results in Table 5.5. Comparison of the results reveals that in the Dreisbach
watershed, BMPs selected and placed by optimization (case 3) would yield nearly three times
better benefit-cost ratio and would cost 2.5 times less than the targeting combination, while
providing the same level of protection against phosphorus and providing even higher
protection against sediment and nitrogen pollution.
Table 5.6 Results for optimization cases
Watershed Variable
Sediment allowable
maximum
o
ctf
'S
Q
_, , allowable
Phosphorus
maximum
,T. allowable
Nitrogen
maximum
Objective function
Watershed cost
„ ,. allowable
Sediment
maximum
PH
^
'e
GO
_, , allowable
Phosphorus
maximum
,T. allowable
Nitrogen
maximum
Objective function
Watershed cost
- From targeting strategy summarized in
- Estimated from SWAT simulations.
Symbol Units
u
ymaxs
ytp
ymaxp
ytn
Z
Ct
yt
, /
ymaxs
ytp
ymaxp
ytn
ymaxn
z
ct
Table 5.5.
t/ha/m
t/ha/m
kg/ha/m
kg/ha/m
kg/ha/m
kg/ha/m
$/$
$
t/ha/m
t/ha/m
kg/ha/m
kg/ha/m
kg/ha/m
kg/ha/m
$/$
$
Optimization
case 3
0.17 -
0.06
0.15 -
0.15
2.10 -
1.55
0.36
165,370
T3
fl}
ta
**r\
UJJ
g
O
c
oj
^
case 4
-
0.085
0.082
1.00
0.21
414,690 -
_
0.48
1.69
26.3
3.12
60,610 -
73
-------
Results for case 3 simulations in the Dreisbach watershed are shown in Figure 5.4(a). In this
case, only maximum monthly sediment and nutrient loads were constrained to their allowable
values obtained from the targeting strategy. A total number of 150 optimization generations
each with a population of 20 strings, were computed. In the top panel in Figure 5.4(a), the
left y-axis reflects benefit-cost ratio for all model evaluations (dots), and right y-axis is total
cost of implementation of the best solution in each optimization generation (dashed line). The
first generation represents a random selection and placement of BMPs, while the last
generation shows the near optimum solution. It is evident that maximum and median fitness
of generations improved as optimization progressed to future generations. This points to the
efficiency of the developed algorithm. Accordingly, the total watershed cost associated with
the best solutions of GA generations generally reduced in successive generations. The fitness
of the best solution of GA generations did not improve once one of the pollutant constituents
reached its allowable value.
w
<3*neratkifl Nwnfcer
50 100
(b)
1 Best individual at Eadi Gsna-alloo
Generation Number
50 109
Cos* (Secondary Axis)
Ge:ier;i;ion Number
SO 100 150
0.1
0 1000 20W 3000
Number of Individual Evaluations
0 1000 200D 3*000
Number of Indr.idual Evaluations
0 1000 MOO 3000 4000 5000
Number of IndMdual Evaluations
-- Sediment (tfliahionth)
Total P (k-gyha/mwrth)
Ml.
— Total N (kgftiaJhionth, Secondary Axis)
1.05 2
,30
0.1
50 100
Gsneratisn Number
50 100 150 209 250
Generation Number
Figure 5.4 (a) Case3 optimization outputs for the Dreisbach watershed over 2000-2009
period, (b) Case 4 optimization outputs for the Dreisbach watershed over 2000-2009 period.
(c) Case 4 optimization outputs for the Smith Fry watershedover2000-2009period.
74
-------
The bottom panel in Figure 5.4(a) demonstrates phosphorus and nitrogen constraints for the
best solution of optimization generations. The values in the figure reflect maximum monthly
pollutant loads simulated over the 2000-2009 period. At the near optimal solution (generation
150), maximum monthly phosphorus load (ymaxp) approached its allowable value (0.15
kg/ha). Estimated maximum monthly nitrogen load (1.55 kg/ha) and maximum monthly
sediment load (0.12 t/ha) were well below allowable values. This indicated that the GA
converged to a near optimal solution with phosphorus being the restricting pollutant. Note
that exceedence of only one of the pollutant loads over its allowable value was sufficient to
penalize a solution. Hence, once monthly phosphorus load approached its allowable value,
regardless of sediment and nitrogen loads being below their allowable values, the fitness of
subsequent optimization generations did not improve, indicating that a near optimal solution
was identified.
Figure 5.4(b) depicts the results of GA computations for case 4, where only the cost
constraint of Eq 5.13 was imposed. The aim of this case was to derive an optimization plan
that would cost the same as the targeting plan. Therefore, the algorithm was designed to
terminate once the total cost of the optimization plan reached the total cost of the targeting
plan ($414,690). A summary of the results for case 4 computations is provided in Table 5.6.
It appeared that for the same watershed cost, the solution from optimization would yield
nearly half the maximum monthly sediment and nutrient loads. Also, the benefit-cost ratio
for near optimal solution obtained from case 4 computations was nearly two times higher
than the one corresponding to the targeting solution.
In the Smith Fry watershed, implementation of BMPs allocated by targeting as shown in
Figure 2.1 would cost $60,610, while reducing average monthly sediment, phosphorus, and
nitrogen loads by nearly 10%, 4.5%, and 3.9%, respectively (see Table 5.5). Since reduction
of pollutant loads was not significant in comparison to the Dreisbach watershed, case 3 was
not investigated for the Smith Fry watershed. The optimization tool was applied to
investigate only case 4, for which only cost constraint was imposed. The analysis aimed at
identification of near optimal combination of BMPs that would cost no more than $60,610,
equal to the cost of implementation of the targeting plan.
75
-------
Figure 5.4(c) shows the results of optimization evaluations for case 4 in the Smith Fry
watershed. Results, also summarized in Table 5.6, indicate that the estimated reduction rate
(Eq 5.18) of maximum monthly sediment, phosphorus, and nitrogen loads from the
optimization procedure would be nearly 5 times higher than the ones estimated for the
targeting plan. Similarly, the benefit-cost ratio would be more than two fold.
A sample of spatial allocation of the optimally placed BMPs in the study watersheds is
depicted in Figure 5.5. In Dreisbach, case 4 was more expensive than case 3 with stricter
sediment and nutrient constraints that were met by allocating additional field borders and
also a few parallel terraces. Field borders appeared to be more widely allocated in the
optimization plans. This is perhaps the case because the cost of implementation of field
borders is significantly less than other types of BMPs studied here (see Table 5.1).
Conversely, grade stabilization structures that are relatively more expensive than the others
were not prescribed in any of the plans. Interestingly, a grassed waterway was delineated at
the very downstream segment of the channel network in all three optimization plans
performed for the study watersheds.
(c)
i Field Boundary
Streams
Grassed Waterway
Parallel Terrace
rq Field Border
1,000
2,000
-1,000
iMeiers
Figure 5.5 Spatial allocation of BMPs from: (a) case 3 computations for Dreisbach
watershed; (b) case 4 computations for Dreisbach; (c) case 4 computations for Smith Fry
watershed.
76
-------
5.6. Discussion
The developed optimization model shows promise for developing watershed restoration and
management plans. This watershed-scale model is well suited to establish relationships
between agricultural practices (particularly structural BMPs), watershed health, and
satisfying aggregate policy objectives. Figure 5.6 shows a sample of reduction of pollutant
loads estimated for the best solution of optimization generations versus their associated
watershed cost for the Dreisbach watershed. Reductions of sediment, phosphorus, and
nitrogen loads at the outlet were weighted by applying weighting factors in Table 5.2 and
were rescaled to [0,1] range such that the maximum value was 1.0. These results clearly
illustrate the capacity of the optimization model to handle the tradeoff amongst
environmental and economic criteria in the objective function. Generating several maps
similar to those of Figure 5.5 enhances the capacity of local managers and any organization
undertaking a watershed planning effort to evaluate attainability of designated water quality
targets at various costs. The foregoing results indicate that the optimization model is likely to
converge to near optimal watershed plans that would achieve the same level of protection
against NFS as targeting strategies at significantly lower costs.
The optimization model benefits from the ability to incorporate spatial heterogeneity in
evaluation of both cost and benefit of BMPs. Critical source areas where BMPs yield utmost
1
o
'= 09
0)
cr
| 0.8
CD
1 0.7
o
Q.
-s 0.6
0
v>
m
ce
n r
•
. . •
.
.
.. *
... *
* ^progress of
. * optimization
*• procedure
+
.*
200,000
400,000
Cost (I)
600,000
800,000
Figure 5.6 Tradeoff between economic criterion (cost) and environmental criteria (weighted
pollutant load reduction) corresponding to case 3 optimization computations for Dreisbach
watershed. Weighted pollutant loads were rescaled to [0,1] range such that the maximum
value was 1.
77
-------
benefits are implicitly identified as the GA progresses to the next generations. Therefore,
benefits of BMPs implemented in critical areas are determined as a function of landscape
characteristics. Although unit cost of implementation of BMPs is the same throughout the
watershed, the total cost is also a function of topographic, soil and land use properties.
Moreover, the optimization procedure handles the interaction between BMPs. The results of
investigated cases in Figure 5.5 indicated that grassed waterways that are implemented
within the channel network were prescribed in near optimal schemes. Conversely, grade
stabilization structures, while very effective in reducing sediment loads (Arabi et al., 2004),
were not allocated in the watersheds, perhaps because of relatively high unit cost of
establishment (see Table 5.1).
The optimization procedure can be performed on a daily, monthly, or annual basis for as
many pollutant constituents as desired. This spatial optimization tool can be linked with other
NFS models through a similar methodology to simulate for constituents that currently cannot
be simulated by SWAT.
The major concern with regard to practicality of the proposed optimization framework is the
CPU time required for computations. The search for near optimal solutions is
computationally demanding mainly because of recurring NPS model simulations. Although
the experienced gained from this research provides guidance as to how select the GA
operating parameters, future studies shall focus on exploring the means to expedite the
procedure. Parallel computing techniques hold great promise in this regard. Also,
improvements in the structure of NPS model(s) that simulates pollutant loads and BMP
representation method would enhance the accuracy of near optimal management plans.
5.7. Conclusions
A GA-based optimization procedure was developed for selection and placement of BMPs.
The sensitivity of the model to different combinations of GA operating parameters, including
population size and replacement rate, was tested in order to identify the most efficient
combination that converges rapidly for a given runtime. For two small watersheds in Indiana,
a setup with a higher number of generations and lower population size was more efficient.
78
-------
However, these results may be site-specific and vary for watersheds with different spatial
scale and characteristics.
The cost effectiveness of the optimized BMPs was compared to that of BMPs prescribed
through targeting strategies in the study watersheds. Targeting strategies are developed based
on field and modeling studies and intend to allocate BMPs in portions of the watershed that
are deemed to extensively contribute to pollutant loads. In the absence of systematic
approaches, however, identifying such locations and also designing most effective BMP
combinations become infeasible at the watershed scale. It was demonstrated that the BMPs
from optimization would achieve the same level of sediment and nutrient reductions with
nearly one third of the cost required for implementation of the targeting scenario.
Conversely, it was shown that an optimized management scheme would likely provide nearly
twice the protection against sediment and nutrient loads for the same amount of money that
would be spent for implementation of the targeting plans in these watersheds.
In all optimization plans, a grassed waterway was allocated to the very downstream of the
channel network and immediately upstream of the watershed outlet where water quality
standards were imposed. Field borders were the most commonly allocated BMPs, perhaps
because of their relatively lower unit cost. Conversely, grade stabilization structures that cost
significantly more than the other types of BMPs in this study were not included in any of the
optimized plans.
79
-------
SECTION 6. CONCLUSIONS
A computational framework was developed and tested on two small watersheds in Indiana to
help watershed management agencies identify effective management practices and their
spatial allocations for sediment and nutrient control. The method utilized simulations of a
distributed-parameter watershed model, Soil and Water Assessment Tool (SWAT), in
conjunction with statistical and optimization techniques.
In this study, the results of a previous study by Arabi et al. (2004) in the same study
watersheds were adopted. First, the SWAT model was calibrated and validated for the study
watersheds. The results indicated that SWAT streamflow, sediment, total P and total N
simulations adequately represent observations. Two error statistics, coefficient of
determination and Nash-Sutcliff efficiency coefficient, were used for calibration purposes.
Model calibration was performed until a coefficient of determination greater than or equal to
0.6 and a Nash-Sutcliff coefficient greater than or equal to 0.5 was obtained for each
constituent. Second, SWAT model parameters were used to represent BMPs based on their
functionality and hydrologic and water quality processes that are modified by their
implementation. Field borders, parallel terraces, grassed waterways, and grade stabilization
structures were the specific BMPs modeled in this study. It would appear that SWAT model
is an appropriate model for representation of agricultural management scenarios. The
effectiveness of BMPs was evaluated by comparing model simulations for two scenarios
representing conditions with and without the presence of BMPs. It was observed that
implementation of the BMPs in the Dreisbach watershed would result in a significant
reduction of sediment and nutrient yields. These reductions would be mainly due to
implementation of grade stabilization structures and grassed waterways that appreciably
reduce the transport capacity of channel network.
A systematic method based on the Generalized Likelihood Uncertainty Estimation (GLUE)
methodology was developed for the estimation of uncertainty bounds for model predictions.
80
-------
The procedure can be utilized for quantification of the margin of safety (MOS) in the TMDL
allocation formula as an explicit uncertainty analysis method. In this study, the Generalized
Likelihood Uncertainty Estimation method was utilized on the two Indiana watersheds. The
forgoing results indicated that the GLUE methodology is suitable not only for determination
of a numeric value for the MOS, but also for addressing the issues that render utility of the
conventional approach for evaluation of BMPs subjective to the questions of uniqueness of
calibrated parameter set, identifiably, and equifinality as described in Beven and Binely (
1992).
Additionally, a computational framework was developed in which investigation of
uncertainty provides complementary quantitative and qualitative information in support of
management and decision making. The analysis focused specifically on two issues. First, the
univariate Regionalized Sensitivity Analysis (RSA) ranking method was applied in
conjunction with the multivariate Tree Structured Density Estimation (TSDE) analysis to
identify critical hydrologic and water quality processes that control the ability to attain
specified water quality conditions. The implications of such analysis provide benefits for
watershed management programs such as TMDL. Application of the RSA-TSDE procedure
for two small watersheds in Indiana, Dreisbach (6.25 km2) and Smith Fry (7.3 km2), showed
that fluvial channel processes appeared to control sediment yield at the outlet of the
watersheds. Therefore, implementation of within-channel BMPs such as grassed waterways
and grade stabilization structures are likely more effective for reducing sediment loads at the
outlet of the study watersheds. However, total nitrogen yield at the outlets was controlled by
upland nitrogen loading, indicating that implementation of parallel terraces or fertilizer
application strategies would be more effective. Also, EPA has been advised to promote
research on translating effectiveness needs into Use Attainability Analysis (UAA)
implications (NRC, 2001). The use of UAA, which is applied for setting new water quality
standards or revising the existing ones, has been restricted due to inadequate technical
guidance from EPA. The results of the TSDE method can provide a numerical probability
value that can be utilized to answer pragmatic questions such as how feasible a target value is
for a desired constituent.
81
-------
Complementary to the RSA-TSDE approach that facilitates identification of effective
management practices, a GA-based optimization model was developed for BMP placement
within watersheds. The cost effectiveness of the optimized BMPs for the study watersheds
was compared to a combination of BMPs that was implemented nearly 30 years ago through
targeting. Results indicated that the benefit-cost ratio of the combination from optimization
was more than two times higher than the one from targeting, providing the same sediment
and nutrient reduction at the outlet. The optimization procedure can be utilized to produce a
number of near optimal solutions, providing flexibility for management and decision making
to meet water quality target values in an economic manner. Also, the sensitivity of the model
to different combinations of operating parameters including population size and replacement
rate was tested in order to identify the most efficient combination which converges faster for
a given runtime. For two small watersheds in Indiana, it would appear that a setup with a
higher number of generations and lower population size was more efficient.
While the foregoing results may be site-specific and limited to watersheds with spatial scale
and characteristics similar to the watersheds used in this study, the developed sensitivity
analysis, uncertainty estimation, and optimization procedures can be applied to in other
watershed studies.
82
-------
BIBLIOGRAPHY
Arabi, M., Govindaraju, R.S., Hantush, M.M. (2007a). A probabilistic approach for analysis
of uncertainty in evaluation of watershed management practices." Journal of Hydrology,
333:459-471.
Arabi, M., R.S. Govindaraju, B. Engel, Hantush, M.M. (2007b). "Multiobjective sensitivity
analysis of sediment and nutrient processes with a watershed model." Water Resources
Research, 43, W06409, doi:10.1029/2006WR005463.
Arabi, M., Govindaraju, R.S., Hantush, M.M. (2006a). "Cost-effective allocation of
watershed management practices using a genetic algorithm." Water Resources Research,
42, W10429, doi:10.1029/2006WR004931.
Arabi, M., Govindaraju, R.S., Hantush, M.M. (2006b). "Role of watershed subdivision on
evaluation of long-term impact of best management practices on water quality." Journal
of American Water Resources Association, 42(2): 513-528.
Arabi, M., Govindaraju, R.S., Hantush, M.M. (2004). "Source identification of sediments and
nutrients in watersheds-Final report." EPA/600/R-05/080, National Risk Management
Research Laboratory, Office of Research and Development, U.S. Environmental
Protection Agency, Cincinnati, OH 45268, 120 pp. Available online at:
www.epa.gov/ORD/NRMRL/pubs/600r05080/600r05080.pdf
Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R. (1998). "Large area hydrologic
modeling and assessment part I: model development." Journal of the American Water
Resources Association, 34(1): 73-89.
Arnold, J.G. and Fohrer, N. (2005). "SWAT2000: Current capabilities and research
opportunities in applied watershed modeling." HydrologicalProcesses, 19(3):563-572.
83
-------
ASAE, 2003. Design, Layout, Construction and Maintenance of Terrace Systems. ASAE
Standards S268.4 FEB03, St. Joseph, MI.
Bastidas, L.A. (1998). "Parameter estimation for hydrometeorological models using
multicriteria methods." Ph.D. dissertation, Department of Hydrology and Water
Resources, University of Arizona, Tuscan, AZ.
Bastidas, L.A., Gupta, H.V., Sorooshian, S., Shuttleworth, W.J., Yang, Z.L. (1999).
"Sensitivity analysis of a land surface scheme using multicriteria methods." Journal of
Geophysical Research, 104(D 16): 19,481 -19,490.
Beck, M.B. (1987). "Water quality modeling: a review of the analysis of uncertainty." Water
Resources Research, 23(8): 1393-1442.
Bekele, E.G. and Nicklow, J.W. (2005). "Multiobjective management of ecosystem services
by interactive watershed modeling and evolutionary algorithms." Water Resources
Research, 41, W10406, doi:10.1029/2005WR004090.
Benaman, J. and Shoemaker, C.A. (2004). "Methodology for analyzing range of uncertain
model parameters and their impact on total maximum daily load process." Journal of
Environmental Engineering, ASCE, 130(6): 648-656.
Beven, KJ. and Binely, A.M. (1992). "The future of distributed models: model calibration
and uncertainty prediction." Hydrological Processes, 6: 279-298.
Beven, KJ. and Freer, J. (2001). "Equifmality, data assimilation, and uncertainty estimation
in mechanistic modeling of complex environmental systems using the GLUE
methodology." Journal of Hydrology, 249(2001): 11-29.
Bracmort, K.S., Arabi, M., Frankenberger, J.R., Engel, B.A., and Arnold, J.G. (2006).
"Modeling long-term water quality impact of structural BMPs." Transactions of the
ASABE, 49(2): 367-374.
84
-------
Buckner, E.R. (2001). "An evaluation of alternative vegetative filter strip models for use on
agricultural lands of the Upper Wabash river watershed." Ph.D. Dissertation, Purdue
University, West Lafayette, IN 47907.
Cangelosi, A., Wether, R., Taverna, J., and Cicero, P. (2001). Revealing the economic value
of protecting the Great Lakes, Northeast-Midwest Institute and the National Oceanic and
Atmospheric Administration, Washington, DC.
Deb, K. (2001). Multi-Objective Optimization Using Evolutionary Algorithms, Wiley,
Chichester, NY.
Demarty, J., Ottle, C., Braud, I, Olioso, A., Frangi, J.P., Gupta, H.V., Bastidas, L.A. (2005).
"Constraining a physically based soil-vegetation-atmosphere transfer model with surface
water content and thermal infrared brightness temperature measurements using a
multiobjective approach." Water Resources Research, 41, W01011, 15 pp.
Dilks, D.W. and Freedman, P.L. (2004). "Improved consideration of the Margin of Safety in
Total Maximum Daily Load development." Journal of Environmental Engineering,
ASCE, 130(6):690-694.
Edwards, D.R., Daniel, T.C., Scott, H.D., Murdoch, J.F., Habiger, M.J., Burkes, H.M.
(1996). "Stream quality impacts of Best Management Practices in a Northwestern
Arkansas basin." Journal of the American Water Resources Association, 32(3): 499-509.
EPA (2003). Introduction to the Clean Water Act, Available online at:
http://www.epa.gov/watertrain/cwa/.
EPA (2005a). "Chapter 7. Analyze data to characterize the watershed and pollutant sources."
in Draft Handbook for Developing Watershed Plans to Restore and Protect Our Waters,
EPA 841-B-05-005, October 2005, Available online at:
http://www.epa.gov/owow/nps/watershed_handbook/.
EPA (2005b). 2002 Section 303(d) List Fact Sheet for INDIANA, Available online at:
http://oaspub.epa.gov/waters/state_rept.control?p_state=IN, accessed on August 2005.
85
-------
Fang, F. and Easter, W.E. (2003). "Pollution trading to offset new pollutant loadings - A
case study in the Minnesota River basin." American Agricultural Economics Association
Annual Meeting, Montreal, Canada.
Fiener, P. and Auerswald, K. (2006). "Seasonal variation of grassed waterway effectiveness
in reducing runoff and sediment delivery from agricultural watersheds in temperate
Europe." Soil and Tillage Research, 87: 48-58.
Fiener, P. and Auerswald, K. (2003). "Concept and effects of a multi-purpose grassed
waterway." Journal of Environmental Quality., 32:927-936.
Freer, J. and Beven, K. (1996). "Bayesian estimation of uncertainty in runoff prediction and
the value of data: An application of the GLUE approach." Water Resources Research,
32(7): 2161-2173.
Green, W.H. and Ampt, G.A. (1911). "Studies on soil physics: 1. the flow of air and water
through soils." Journal of Agricultural Sciences, 4: 11-24.
Griffin, C.B. (1995). "Uncertainty analysis of BMP effectiveness for controlling nitrogen
from urban nonpoint sources." Journal of the American Water Resources Association,
31(6): 1041-1050.
Heuvelmans, G., Muys, B., Feyen, J. (2004). "Evaluation of hydrological model parameter
transferability for simulating the impact of land use on catchment hydrology." Physics
and Chemistry of the Earth, 29(2002): 739-747.
Ice, G. (2004). "History of innovative best management practice development and its role in
addressing water quality limited waterbodies." Journal of Environmental Engineering,
ASCE, 130(6): 684-689.
Indiana Environmental Quality Incentives Program (2004). 2004 EQIP Cost List Information
Available online at: http://www.in.nrcs.usda.gov/programs/2004eqip/state_costlist.html.
Kirsch, K., Kirsch, A., Arnold, J.G. (2002). "Predicting sediment and phosphorous loads in
the Rock River basin using SWAT." the Transactions ofASAE, 45(6): 1757-1769.
86
-------
Lake, J. and Morrison, J. (1977). "Environmental impacts of land use on water quality, Final
report on the Black Creek project-project data." EPA-905/9-77-007-C, U.S.
Environmental Protection Agency, Chicago, IL.
Lake, J., Morrison, J. (1978). "Environmental impact of land use on water quality- Final
report on the Black Creek project- Technical report." EPA-905/9-77-007-B, U.S.
Environmental Protection Agency, Chicago, IL.
Lenhart, T., Eckhardt, K., Fohrer, N., Frede, H.-G. (2002). "Comparison of two approaches
of sensitivity analysis." Physics and Chemistry of the Earth., 27 (2002): 645-654.
Liu, Y., Gupta, H.V., Sorooshian, S., Bastidas, L.A., Shuttleworth, W.J. (2004). "Exploring
parameter sensitivities of the land surface using a locally coupled land-atmosphere
model." Journal of Geophysical Research, 109, D21101, 13 pp.
McKay, M.D., Beckman, R.J., Conover, W.J. (1979). "A Comparison of Three Methods of
Selecting Values of Input Variables in the Analysis of Output from a Computer Code."
Technometrics, 21: 239-245.
Meixner, T., Gupta, H.V., Bastidas, L.A., Bales, R.C. (1999). "Sensitivity analysis using
mass flux and concentration." HydrologicalProcesses, 13:2233-2244.
Morris, M.D. (1991). "Factorial sampling plans for preliminary computational experiments."
Technometrics, 33: 161-174.
Morrison, J. and Lake, J. (1983). "Environmental impacts of land use on water quality, Black
Creek Project- Final report." Allen County Soil and Water Conservation District, Fort
Wayne, IN.
Mostaghimi, S., Park, S.W., Cook, R.A., Wang, S.Y. (1997). "Assessment of management
alternatives on a small agricultural watershed." Water Research, 31(8): 1867-1878.
Muleta, M.K. and Niklow, J.W. (2005). "Sensitivity and uncertainty analysis coupled with
automatic calibration for a distributed watershed model." Journal of Hydrology, 306
(2005): 127-145.
87
-------
Nash, I.E. and Sutcliffe, J.V. (1970). "River flow forecasting through conceptual models part
1- A discussion of principles." Journal of Hydrology, 10(1970): 282-290.
Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., King ,K.W. (2002). Soil and Water
Assessment Tool Theoretical Documentation, Version 2000. Grassland, Soil and Water
Research Laboratory, Agricultural Research Service, Temple. TX.
NRC (2001). 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.
Osidele, O.O., Zeng, W., Beck, M.B. (2003). "Coping with uncertainty: a case study in
sediment transport and nutrient load analysis." Journal of Water Resources Planning and
Management, 129(4): 345-355.
Pitman, AJ. (1994). "Assessing the sensitivity of land-surface scheme to the parameter
values using a single column model." Journal of Climate, 7(12)1856-1869.
Ribaudo, M.O., Colacicco, D., Barbarika, A., and Young, C.E. (1989). "The economic
efficiency of voluntary soil conservation programs." Journal of Soil and Water
Conservation, 44: 40-43.
Ringuest, J.L. (1992). Multiobjective Optimization: Behavioral and Computational
Considerations, Boston, MA: Kluwer.
Saleh, A., Arnold, J.G., Gassman, P.W., Hauch, L.M., Rosenthal, W.D., Williams, J.R.,
McFarland, A.M.S. (2000). "Application of SWAT for the Upper North Bosque River
watershed." the Transactions of ASAE, 43(5): 1077-1087.
Saltelli, A., Chan, K., Scott, E.M. (2000). Sensitivity Analysis, John Wiley & Sons Ltd., West
Sussex PO19 8SQ, England.
-------
Santhi, C., Arnold, J.G., Williams, J.R., Dugas, W.A., Srinivasan, R., Hauck, L.M. (2001).
"Validation of the SWAT model on a large river basin with point and nonpoint sources."
Journal of the American Water Resources Association., 37(5): 1169-1188.
Santhi, C., Srinivasan, R., Arnold, J.G., Williams, J.R. (2003). "A modeling approach to
evaluate the impacts of water quality management plans implemented in the Big Cypress
Creek Watershed." second conference on Watershed Management to Meet Emerging
TMDL Environmental Regulations, Albuquerque, NM, pp. 384-394.
Spear, R.C. and Hornberger, G.M. (1980). "Eutrophication in Peel Inlet II: Identification of
critical uncertainties via generalized sensitivity analysis." Water Research, 14: 43-49.
Spear, R.C., Grieb, T.M., Shang, N. (1994). "Parameter uncertainty and interaction in
complex environmental models." Water Resources Research, 30(11): 3159-3169.
Srivastava, P., Robillard, P.O., Hamlet, J.M., and Day, R.L. (2002). "Watershed optimization
of best management practices using AnnAGNPS and a genetic algorithm." Water
Resources Research, 38(3): 1021.
USDA Soil Conservation Service (1972). "Section 4- Hydrology, Chapters 4-10." National
Engineering Handbook, USDA, Washington, DC.
USDA (2003). The 2002 Farm Bill: provisions and economic implications, Available online
at: http://www.ers.usda.gov/Features/FarmBill/.
USDA-NRCS (2005). "National Conservation Practice Standards- NHCP." Electronic Field
Office technical Guide (eFOTG), Section IV-Practice Standards and Specifications,
Available online at: http://www.nrcs.usda.gov/Technical/Standards/nhcp.html, Accessed
on August 2006.
van Griensven, A. and Bauwens W. (2003). "Multiobjective autocalibration for
semidistributed water quality models." Water Resource Research, 39(12): 1348,
doi: 10.1029/2003WR002284, 2003.
89
-------
van Griensven, A., Meixner, T., Grunwald, S., Bishop, T., Diluzio, M., Srinivasan, R. (2006).
"A global sensitivity analysis tool for the parameters of multi-variable catchment
models." Journal of Hydrology, 324(2006): 10-23.
Vazquez-Amabile, G., Engel, B.A., Flanagan, D.C. (2006). "Modeling and risk analysis of
nonpoint-source pollution caused by atrazine using SWAT." Transactions of the ASABE,
49(3): 667-678.
Veith, T.L., Wolfe, M.L., and Heatwole, C.D. (2004). "Cost-effective BMP placement:
optimization versus targeting." Transactions of the ASAE, 47(5): 1585-1594.
White, K.L. and Chaubey, I. (2005). "Sensitivity analysis, calibrations, and validation for a
multisite and multivariate SWAT model." Journal of the American Water Resources
Association, 41 (5): 1077-1089.
Whittaker, G., Fare, R., Srinivasan, R. and Scott, D.W. (2003). "Spatial evaluation of
alternative nonpoint nutrient regulatory instruments." Water Resources Research, 39(4)
1079, doi:l0.1029/2001 WR001119, 2003.
Williams, J.R. (1969). "Flood routing with variable travel time or variable storage
coefficients." the Transactions of the ASAE, 12(1): 100-103.
Williams, J.R. (1975). "Sediment-yield prediction with universal equation using runoff
energy factor, in Present and prospective technology for predicting sediment yield and
sources." Proceedings of the Sediment Yield Workshop, Oxford, MS, pp. 244-252.
Winston, W.L. and Venkataramanan, M. (2003). Introduction to Mathematical Programming
Brooks/Cole, California 93950, USA.
Wischmeier, W.H. and Smith, D.D. (1978). "Predicting Rainfall Erosion Losses, a Guide to
Conservation Planning." Agriculture Handbook NO. 537, US Department of Agriculture,
Washington D.C.
90
-------
Zhang, H.X. and Yu, S.L. (2004). "Applying the first-order error analysis in determining the
margin of safety for total maximum daily load computations." Journal of Environmental
Engineering, 130(6): 664-673.
Zheng, Y. and Keller, A.A. (2006). "Understanding parameter sensitivity and its
management implications in watershed-scale water quality modeling." Water Resources
Research, 42, W05402, 14 pp.
Zitzler, E. and Thiele, L. (1999). "Multiobjective evolutionary algorithms: a comparative
case study and the Strength Pareto approach." IEEE Transactions on Evolutionary
Algorithms, 3(4): 257-271.
91
-------
SEPA
United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Cincinnati, OH 45268
Official Business
Penalty for Private Use
$300
Contract Number: 4C-R330-NAEX
January 2007
------- |