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 •".'- ' J *5'3 T=\ Ay = y(x,t,td,a) Eq 5.4 where AyiT is annual reduction of constituent /' (i.e. sediment in t/yr, phosphorus in kg/yr, and nitrogen in kg/yr) for yearr at the outlet of the study watershed, td (yr) is BMP design life, and bm is the monetary benefit of BMP as a result of reducing constituent /' from Table 5.2. In 5.3, x is the state of the watershed at any given time t,a denotes the watershed management plan, and function y reflects mathematical relationships in the SWAT model that are used for representation of hydrologic and water quality processes. 63 ------- 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 ------- |