Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0047356/00001
## Material Information- Title:
- Hydrologic and Water Quality Model Reliability With Global Sensitivity Analysis Improvements and Applications
- Creator:
- Khare, Yogesh P
- Place of Publication:
- [Gainesville, Fla.]
Florida - Publisher:
- University of Florida
- Publication Date:
- 2014
- Language:
- english
- Physical Description:
- 1 online resource (209 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Agricultural and Biological Engineering
- Committee Chair:
- MARTINEZ,CHRISTOPHER J
- Committee Co-Chair:
- MUNOZ-CARPENA,RAFAEL
- Committee Members:
- JONES,JAMES W
KAPLAN,DAVID A BOTTCHER,ADELBERT B - Graduation Date:
- 12/19/2014
## Subjects- Subjects / Keywords:
- Arid zones ( jstor )
Groundwater ( jstor ) Hydrological modeling ( jstor ) Mathematical independent variables ( jstor ) Modeling ( jstor ) Parametric models ( jstor ) Sensitivity analysis ( jstor ) Soils ( jstor ) Trajectories ( jstor ) Watersheds ( jstor ) Agricultural and Biological Engineering -- Dissertations, Academic -- UF hydrology -- model -- reliability -- sensitivity -- water-quality Lake Okeechobee ( local ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) government publication (state, provincial, terriorial, dependent) ( marcgt ) born-digital ( sobekcm ) Electronic Thesis or Dissertation Agricultural and Biological Engineering thesis, Ph.D.
## Notes- Abstract:
- Hydrologic and water quality (H/WQ) models can well-assess the impacts of natural and anthropogenic disturbances on natural systems. Today, these models have become extremely important policy making tools for tackling a range of water resources problems e.g. Best Management Practices, Total Maximum Daily Loads, etc. For well-informed decision making, applications of these models should undergo thorough reliability/global sensitivity and uncertainty analyses (GSA/UA) evaluation. Though the H/WQ modeling community has begun to realize the importance of GSA/UA, they are often ignored in practice. This study aimed to understand and improve SA/UA techniques and their applications to H/WQ and drought models. This work is divided into three specific research objectives. A drought index, the Agricultural Reference Index of Drought (ARID), a model with correlated parameters was used as a test case for the first objective. An innovative approach was proposed to incorporate parameter correlation while performing GSA using a variance based technique. GSA results obtained using this approach were compared to those of a correlation-based SA methodology. Results were comparable. ARID was found to be monotonic and additive in nature and was most sensitive to the root zone depth and soil hydraulic properties. The second part of this work proposed a new multi-criteria trajectory-based sampling strategy for the screening SA method of Elementary Effects (EEs). The concept of generating uniformly distributed parameter samples was merged with the often-used goal of maximizing sample spread. The evaluation of this strategy with benchmark strategies indicated that it was more than an order of magnitude computationally faster for high dimensional models and that it also somewhat improves parameter screening. The last objective was aimed at the evaluation of the routing and pollutant attenuation module of the Watershed Assessment Model (WAM). GSA/UA showed that WAM was most sensitive to empirical/conceptual parameters, indicating their importance in the model calibration/validation process and the need for monitoring efforts to verify the selection of these parameters. As the first formal evaluation of WAM by GSA/UA, the results obtained in this work are a valuable contribution to the application of WAM to addressing water quality issues in the state of Florida. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2014.
- Local:
- Adviser: MARTINEZ,CHRISTOPHER J.
- Local:
- Co-adviser: MUNOZ-CARPENA,RAFAEL.
- Electronic Access:
- RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2015-12-31
- Statement of Responsibility:
- by Yogesh P Khare.
## Record Information- Source Institution:
- UFRGP
- Rights Management:
- Applicable rights reserved.
- Embargo Date:
- 12/31/2015
- Resource Identifier:
- 974372570 ( OCLC )
- Classification:
- LD1780 2014 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

PAGE 1 HYDROLOGIC AND WATER QUALITY MODEL EVALUATION WITH GLOBAL SENSITIVITY ANALYSIS : IMPROVEMENTS AND APPLICATIONS By YOGESH P . KHARE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2014 PAGE 2 Â© 2014 Yogesh P . Khare PAGE 3 To my m om PAGE 4 4 ACKNOWLEDGMENTS I would like to thank my advisor Dr. Christopher Martinez for his constant support and encouragement over the past five years. I could not have achieved this goal without Carpena and all the members of the graduate committee: Drs Jim Jones, Del Bottcher, and David Kaplan. I would also like to thank Dr. Andrew James from the Soil and Water Engineering Technology Inc. (SWET) for his help in understanding the Watershed Assessment Model (WAM), the University of Florida (UF) Research Computing Center team, and all funding agencies : the UF Center for Landscape Conservation and Ecology (CLCE), the Florida Sea Grant (FSG) and others, for their support in making this research possible. Special thanks to Dr. Ashish Mehta and his family for their love and guidance since I began my journe y of graduate studies at UF. I also thank my friends and well wishers in India, in Gainesville, and in the Agricultural and Biological Engineering Department for their direct or indirect support to make this journey smooth. I am very grateful to my parents , Prakash and Sunita Khare, my sisters, Vrushali and Seema, and my in laws, Uday and Smita Chitale, for their love, patience, and understanding my educational priorities. Last, but certainly not least, I thank my wife, Janhavi, for being my strength, motiv ation, and the better half of me. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ ............... 4 LIST OF TABLES ................................ ................................ ................................ ........................... 8 LIST OF FIGURES ................................ ................................ ................................ ....................... 10 ABSTRACT ................................ ................................ ................................ ................................ ... 13 CHAPTER 1 INTRODUCTIO N ................................ ................................ ................................ .................. 15 Research Trends in Hydrologic and Water Quality (H/WQ) Modeling ................................ . 15 Sensitivity and Uncertainty Analysis: An Overview ................................ .............................. 18 Sensitivity Analysis Techniques ................................ ................................ ...................... 18 One at a time ................................ ................................ ................................ ............ 19 Method of e lementary effects/ Morris method ................................ ......................... 19 Correlation/regression methods ................................ ................................ ................ 21 Regionalized sensitivity analysis (RSA)/ Hornberger Spear Young method .......... 22 Multi objective regionalized sensitivity analysis (MORSA) ................................ ... 23 Variance based methods ................................ ................................ ........................... 23 Recent advancements in SA techniques ................................ ................................ ... 26 Uncertainty Analysis Techniques ................................ ................................ .................... 26 Sensitivity and Uncertainty Analysis in H/WQ Modeling: Applications and Issues ............. 28 Objectives ................................ ................................ ................................ ............................... 31 2 PARAMETER VARIABILITY AND DROUGHT MODELS A STUDY USING THE AGRICULTURAL INDEX FOR DROUGHT (ARID) ................................ ................ 41 Materials and Methods ................................ ................................ ................................ ........... 44 Model ................................ ................................ ................................ ............................... 44 Reference evapotranspiration (ETo) ................................ ................................ ........ 45 Evapotranspiration (ET) ................................ ................................ ........................... 46 Surface runoff (R) ................................ ................................ ................................ .... 46 Deep drainage (D) ................................ ................................ ................................ .... 47 Application/ Case Study Details ................................ ................................ ...................... 47 GSA/UA Procedure ................................ ................................ ................................ ......... 47 Simulation Set Up ................................ ................................ ................................ ........... 48 Input data and parameter PDFs ................................ ................................ ................ 48 GSA/UA Methods ................................ ................................ ................................ .... 51 Model Output, Uncertainty and Sensitivity Indicators ................................ ............. 52 Results ................................ ................................ ................................ ................................ ..... 53 Histograms and CDFs of HT (based on model run ....................... 53 Scatterplots (based on model runs for correlation analysis) ................................ ............ 54 PAGE 6 6 Sensitivity Indices ................................ ................................ ................................ ........... 54 Correlation a nalysis based sensitivity indices ................................ .......................... 54 ................................ ................................ ................................ ........... 55 Effect of Location ................................ ................................ ................................ ............ 56 Discussion ................................ ................................ ................................ ............................... 57 Comparison With Other Studies ................................ ................................ ...................... 58 Summary ................................ ................................ ................................ ................................ . 61 3 A MULTI CRITERIA TRAJECTORY BASED PARAMETER S AMPLING STRATEGY FOR THE SCR EENING METHOD OF ELE MENTARY EFFECTS ............ 77 Material and Methods ................................ ................................ ................................ ............. 81 Method of EEs ................................ ................................ ................................ ................. 81 Trajectories Based on Sampling for Uniformity (SU) ................................ .................... 82 Num erical Experiments ................................ ................................ ................................ ... 83 Testing for uniformity of the generated distributions ................................ .............. 84 Time comparison ................................ ................................ ................................ ...... 85 Euclidean distance ................................ ................................ ................................ .... 85 Parameter screening efficiency ................................ ................................ ................ 86 Results ................................ ................................ ................................ ................................ ..... 89 Uniformity of Generated Parameter Distributions ................................ .......................... 89 Sampling Time Requirement ................................ ................................ ........................... 89 Spread of Trajectories/ Euclidean Distances ................................ ................................ ... 90 Parameter Screening Efficiency ................................ ................................ ...................... 90 Discussion ................................ ................................ ................................ ............................... 91 Theoretical Basis: Sampling for Uniformity ................................ ................................ ... 91 Sample Generation Runtime Efficiency ................................ ................................ .......... 92 Spread of Traject ories in Hyperspace and Parameter Screening Performance ............... 93 Other Considerations ................................ ................................ ................................ ....... 95 Summary ................................ ................................ ................................ ................................ . 95 4 GLOBAL SENSITIVITY AND UNCERTAINTY ANALYSIS OF THE WATERSHED ASSESSMENT MODEL (WAM) ................................ .............................. 108 Material and Methods ................................ ................................ ................................ ........... 111 Model ................................ ................................ ................................ ............................. 111 Source cell simulations ................................ ................................ ........................... 112 Flow and nutrient routing ................................ ................................ ....................... 112 Study Watershed and Model Set up ................................ ................................ .............. 113 Evaluation Methodolog ies ................................ ................................ ............................. 115 Global sensitivity analysis ................................ ................................ ...................... 115 Uncertainty analysis and parameter estimation ................................ ...................... 1 17 Experiment Details ................................ ................................ ................................ ........ 118 Parameter distributions ................................ ................................ ........................... 118 Outputs of interest ................................ ................................ ................................ .. 119 Other details ................................ ................................ ................................ ........... 121 Results ................................ ................................ ................................ ................................ ... 121 PAGE 7 7 Global Sensitivity Analysis ................................ ................................ ........................... 121 Elementary Effect sensiti vity analysis ................................ ................................ ... 121 ................................ ................................ ................................ ....... 123 Uncertainty Analysis ................................ ................................ ................................ ..... 124 Monte Carlo Filtering ................................ ................................ ................................ .... 125 Discussion ................................ ................................ ................................ ............................. 126 Implications for WAM Applications, Potential Model and Documentation Improvements ................................ ................................ ................................ ............ 126 Nutrient attenuation ................................ ................................ ................................ 126 Model performance ................................ ................................ ................................ 128 NSE daily flow and D_RO_kDF ................................ ................................ ............ 130 Inconsistencies in WAM documentation ................................ ............................... 131 Evaluation Techniques ................................ ................................ ................................ .. 131 ................................ ................................ .... 131 Misclassification of parameters ................................ ................................ .............. 132 Identification of parameter interaction and non monotonic effects ....................... 133 Normalization of sensitivity measures ................................ ................................ ... 134 Summary ................................ ................................ ................................ ............................... 134 5 CONCLUSIONS ................................ ................................ ................................ .................. 165 APPENDIX A ALGORITHM AND MATLAB PROGRAM FOR THE SAMP LING FOR UNIFORMITY SAMPLING STRATEGY ................................ ................................ .......... 170 B TEST FUNCTIONS ................................ ................................ ................................ ............. 173 C WAM/BLASROUTE PARAMETERS ................................ ................................ ................ 176 LIST OF REFERENCES ................................ ................................ ................................ ............. 187 BIOGRAPHICAL SKETCH ................................ ................................ ................................ ....... 209 PAGE 8 8 LIST OF TABLES Table page 1 1 Summary of widely used watershed scale models ................................ ............................. 33 1 2 Summary of selected sensitivity analysis applications in H/WQ modeling studies .......... 34 1 3 Summary of parametric uncertainty analysis applications in selected H/WQ modeling studies ................................ ................................ ................................ ................ 36 1 4 Summary of input data uncertainty analysis applications in selected H/WQ modeling studies ................................ ................................ ................................ ................................ 37 1 5 Summary of applications of sensitivity analysis techniques in H/WQ modeling pre and post 2000 based on Table 1 2 and Yuan et al. (2014) ................................ ................. 38 2 1 List of input parameters and summary of their estimated/ proposed probability distribution functions ................................ ................................ ................................ ......... 63 2 2 Characte ristics of the distributions of high threat (HT) 1 (LL=lower limit, UL=upper limit) for the four soils used in the study (SL= Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) ................................ ................................ ........ 64 2 3 Partial Correlation Coefficient (PCC), Partial Rank Correlation Coefficient (PRC) and Parameter Rankings ................................ ................................ ................................ .... 65 2 4 ................................ ................................ .............. 66 2 5 Coefficient of model determination values (R 2 ) based on raw and ranked output obtained in partial correlation analysis (for the four soils used in the study (SL = Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) ........... 67 2 6 Summary of wilting point ( wp ), field capacity ( fc ), and Total available water content ( TAWC ) (SL= Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) ................................ ................................ ................................ .......................... 68 2 7 Summary of uncertainty and sensitivity analyses of ARID and other soil w ater balance models ................................ ................................ ................................ ................... 69 3 1 Elementary effect methods and sampling strategies for parameter screening sensitivity analysis ................................ ................................ ................................ ............. 97 3 2 Test function configurations and the number of important parameters (based on analytical expression/ previously calculated values of sensitivity indices) ....................... 98 PAGE 9 9 3 3 The average fraction of parameter distributions generated by the six sampling strategies failing t he Chi square test for discrete uniform distribution at three significance levels for six sampling strategies ................................ ................................ ... 99 3 4 Summary of s creening efficiency skill scores for six sampling strategies: Morris Method, Optimized Trajectories, Modified Optimized Trajectories, and Sampling for Uniformity in three configurations (SU1000, SU300, and SU1) ................................ ..... 100 3 5 Qualitative assessment of sampling strategies used in this study across all evaluation criteria ................................ ................................ ................................ .............................. 102 4 1 Percentage land use distribution for the S 191 watershed ................................ ............... 136 4 2 Uncertainty analysis, and Monte Carlo Filtering ................................ ............................. 137 4 3 Criteria used for Monte Carlo Filtering (MCF) ................................ ............................... 138 4 4 Summary of important WAM/BLASROUTE parameters for flow, BOD, and TSS related outputs. Numbers in parenthesis indicate parameter rank based on * . Please refer Table C 1 for parameter definitions. ................................ ................................ ....... 139 4 5 Summary of important WAM/BLASROUTE parameters for TP and TN related outputs. Numbers in parenthesis indicate parameter rank based on * . Please re fer Table C 1 for parameter definitions. ................................ ................................ ................ 140 4 6 Summary of WAM/BLASROUTE outputs posterior to Monte Carlo Filtering (MCF) . 141 4 7 Summary of parameter ranges pre and post Monte Carlo Filtering (MCF) .................... 142 4 8 Correlation coefficients between ranks of parameters obtained from EE screening ................................ ................................ ................................ ......... 143 4 9 proposed criteria for identification of non monotonic/interaction effect from EE analysis ................................ ................................ ................................ ............................. 144 C 1 Summary of WAM/BLASROUTE parameters and their probabilistic distributions used for the purpose of parameter screening using the method of Elementary Effects ... 176 PAGE 10 10 LIST OF FIGURES Figure page 1 1 Classification of sensitivity analysis techniques ................................ ................................ 39 1 2 Relative usage of various sensitivity analysis techniques applied in H/WQ studies. Results are based on 115 studies ................................ ................................ ........................ 40 2 1 Schematic Soil Water Balance. The square represents the maximum volume of water stored in the soil, with three moisture regions from easily extractable to more difficult (saturation, field capacity, wilting point) ................................ ............................. 70 2 2 Locations of weather stations used in this study ................................ ................................ 71 2 3 Long term mean monthly precipitation, maximum temperature, and minimum temperature at the five study locations (based on 1980 2003 daily data) .......................... 72 2 4 ................................ .............. 73 2 5 Scatterplots of HT against model parameters for the four soil types at five study locations . ................................ ................................ ................................ ............................ 75 2 6 Scatterplots of HT against Total Available Water Content ( TAWC ) for all soil types and location combinations ................................ ................................ ................................ . 76 3 1 Schematic representation of parameter levels and constant offset coordinates/ jump size for an arbitrary parameter ................................ ................................ ......................... 103 3 2 Schematic representation of sample generation using Sampling for Uniformity (SU) for the case of 2 parameters, 4 trajectories . ................................ ................................ ..... 104 3 3 Time required for sample generation of size k at two trajectory configurations ( r = 6 and r = 10) using six parameter samp ling strategies: MM ( Q = r ), OT ( Q = 1000 ), MOT ( Q = 1000), SU1000 ( Q = 1000), SU300 ( Q = 300), a nd SU1 ( Q = 1) ................. 105 3 4 Box plots of the Euclidean distances (measure of tra jectory spread): MM, OT, MOT, SU1000 , SU300 , and SU1 ; across range of model sizes (i.e. number of parameters, k ) for r = 10 . ................................ ................................ ................................ ......................... 106 3 5 Parameter screening efficiencies for six sampling strategies: MM , OT, MOT, SU1000, SU300 , and SU1 ; across all test functions using Approach 1 (a, b) and A pproach 2 (c, d) ................................ ................................ ................................ ............. 107 4 1 Determination of unique model cells in WAM (adopted from WAM Documentation, SWET 2014). ................................ ................................ ................................ ................... 145 PAGE 11 11 4 2 Schematic representation the BUCSHELL module (adopted from WAM Documentation, SWET 2014). ................................ ................................ ......................... 146 4 3 Schematic of the BLASROUTE sub module flow directions and distance grids (adopted from Bottcher et al. 2012). ................................ ................................ ................ 147 4 4 Maps showing (a) Location and topography of S 191 watershed, (b) Land use within S 191 watershed and (c) Distribution of soils within S 191 watershed. .......................... 148 4 5 Elementary Effects sensitivity analysis results for the flow relate d outputs in and spaces . ................................ ................................ ................................ .............. 149 4 6 Elementary Effects sensitivity anal ysis results for TP outputs and spaces . ................................ ................................ ................................ .............................. 150 4 7 Elementary Effects sensitivity anal ysis results for TN outputs spaces. ................................ ................................ ................................ .............................. 151 4 8 Elementary Effects sensitivity analysis results for the BOD related outputs ................................ ................................ ................................ .............. 152 4 9 Elementary Effects sensitivity analysis results for the TSS related outputs ................................ ................................ ................................ .............. 153 4 10 S and ST ), for flow outputs. Parameters for which both S and ST were <0.5% or 0.005 were not plotted. ................................ ................................ ................................ ............................. 154 4 11 S and ST ), for TP outputs. Parameters for whi ch both S and ST were <0.5% or 0.005 were not plotted. ................................ ................................ ................................ ............................. 155 4 12 sensitivity indices ( S and ST ), for TN outputs. Parameters for which both S and ST were <0.5% or 0.005 were not plotted. ................................ ................................ ................................ ............................. 156 4 13 Probability Density Function (PDF), Cumulative Density Function (CDF), 95% confidence interval (2.5% 97.5%), and median values for flow outputs . . ....................... 157 4 14 Pr obability Density Function , C umulative Density Function , 95% confidence interval (2.5% 97.5%), and median values for TP outputs . ................................ ............. 158 4 15 Pr obability Density Function , C umulative Density Function , 95% confidence interval (2.5% 97.5%), and median values for TN outputs. ................................ ............ 159 4 16 Posterior distributions of WAM/BLASROUTE outputs analyzed by M onte C arlo F iltering . ................................ ................................ ................................ ........................... 160 PAGE 12 12 4 17 sensitivity analysis and uncertainty analysis part 1 . ................................ ........................ 161 4 18 sensitivity analysis and uncertainty analysis part 2 ................................ ......................... 162 4 19 sensitivity analysis and uncertainty analysis part 3 ................................ ......................... 163 4 20 Scatterplot of NSE daily flow and D_RO_kDF ................................ ............................... 164 A 1 Flowchart for Sampling for Uniformity (SU) to generate sample ( r trajectories) for k parameters with the resampling size of Q ................................ ................................ ........ 172 C 1 Runoff unit hydrograph options used for the EE analysis. RO_Hydrograph 3 is the default runoff hydrogrpah ................................ ................................ ................................ 184 C 2 Sub surface unit hydrograph options used for the EE analysis. Sub surface Unit Hydrograph 1 is the default groundwater unit hydrograph in WAM. Note that a ll hydrogrpahs are of 90 day duration. For brevity only first 10 days are shown. .............. 185 C 3 Runoff unit hydrograph options used for ................................ ........ 186 PAGE 13 13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy HYDROLOGIC AND WATER QUALITY M ODEL EVALUATION WITH GLOBAL SENSITIVITY ANALYSIS: IMPROVEMENTS AND APPLICATIONS By Yogesh P. Khare December 2014 Chair: Christop her Martinez Cochair: oz Carpena Major: Agricultural and Biological Engineering Hydrologic and water quality (H/WQ) models can well assess the impacts of natural and anthropogenic disturbances on natural systems. Today , these models have become extremely important policy making tools for tackling a range of water resources problems e.g. Best Management Practices, Total Maximum Daily Loads, etc. For well informed decision making, a pplications of these models should undergo th o rough reliability/ gl obal sensitivity and uncertainty analyses (GSA/UA) evaluation . Though the H/WQ modeling community has begun to realize the importance of GSA/UA , they are often ignored in practice. This study aim ed to understand and improv e SA/UA techniques and their applications to H/WQ and drought models. This w ork is divided into three specific research objectives. A drought index , the Agricultural Reference Index of Drought ( ARID ), a model with correlated parameters was used as a test case for the first o bjective. An innovative approach was proposed to incorporate parameter correlation while performing GSA using a variance based technique. GSA results obtained using this approach were compared to those of a correlation based SA methodology. Results were co mparable. ARID was found to be monotonic and additive in nature and was most sensitive to the root zone depth and soil hydraulic properties. PAGE 14 14 T he second part of this work proposed a new multi criteria trajectory based sampling strategy for the screening SA method of Elementary Effects (EEs). The c oncept of generating uniformly distributed parameter samples was merged with the often used goal of maximizing sample spread. The evaluation of this strategy with benchmark strategies indicated that it was more tha n an order of magnitude computationally faster for high dimensional models and that it also somewhat improves parameter screening. The last objective was aimed at the evaluation of the routing and pollutant attenuation module of the Watershed Assessment Mo del (WAM) . GSA/UA showed that WAM was most sensitive to empirical /conceptual parameters , indicating their importance in the model calibration/validation process and the need for monitoring efforts to verify the selection of these parameters . As the first f ormal evaluation of WAM by GSA/UA, the results obtained in this work are a valuable contribution to the application of WAM to addressing water quality issues in the state of Florida . PAGE 15 15 CHAPTER 1 INTRODUCTION Hydrologic and water quality (H/WQ) simulation models replicate a number of natural processes such as surface and sub surface flows, sediment transport, chemical and nutrient flows etc. These process specific models are often integrated together to build c omplex computer programs which can simulate large natural systems such as watersheds. Such large scale models, commonly known as watershed models, can very well assess the i mpacts of natural ( e.g. climate change) and anthropogenic (land use change, defores tation etc.) disturbances on natural systems ( Migliaccio and Srivastava, 2007; Singh and Frevert, 2006 ). Due to this, today these H/WQ models (individual process models as well as watershed models) have become extremely important tools for tackling a range of water resources, environmental , and social problems. Though a number of H/WQ models are available they vary in their abilities due to their inherent properties e.g. scale, underlying fundamental processes, numerical schemes, data requirements, algorithms , purpose etc. Comprehensive summaries of some of the widely used large scale H/WQ mo dels and their properties can be found in Singh and Frevert (2002a, b; 2006), Bora h and Bera (2003, 2004), Daniel et al. (2011) , Sing h (2012) etc. A few popular watershed models are summarized in Table 1 1 including the Watershed Assessment Model WAM (Bottcher et al., 2002 ; Bottcher et al., 2012 ) , the watershed model of interest for a part of this work. Research Trends in Hydrologi c and Water Quality (H/WQ) Mode ling Research foci in H/WQ modeling have shifted over the last five decades. While earlier efforts were focused towards model building, algorithms etc., later, a great amount of attention was paid towards measurement of hydrologic data and developments in data measurement technologies (e.g. remot e sensing and radar technology) and data managemen t systems (e.g. PAGE 16 16 Geographic Information System (GIS) and Database Management System (DBMS) ). With these developments H/WQ models have become very comprehensive but equally complex. Some of the recent research foci include data driven modeling technologies l ike artificial neural networks, fuzzy logic , and genetic algorithms (e.g. Govindaraju and Rao, 2000; Srivastava et al. 2002 etc.). Though these techniques have shown to perform efficiently , there exists an inherent disadvantage/ limitation when the availab le data is scarce or unreliable or totally unavailable (e.g. modeling ungauged watersheds). This , in itself , has attracted considerable attention since the late 1990s (e.g. Wagener et al., 2004; Vogel, 2006) and remains an active area of research. Difficu lties in model calibration and validation arising from the complexities of watershed models have led to developments in optimization algorithms, objective functions, automated parameter estimation and calibration etc. (Sorooshian and Gupta, 1995; Gupta et al., 1999 etc.). However, another group of hydrologists and scientists have realized the importance of the stochastic behavior of model parameters; non unique parameter sets for equivalent model performance (i.e. equifinality) and hence the need of probabi listic approaches while modeling any kind of natural system. An important limitation for this kind of modeling is the intensive need for computational resources. In spite of the above mentioned issues, model calibration and validation has largely remained address this issue the American Society of Agricultural and Biological Engineers (ASABE) has formed a task committee which is working on setting up guidelines for model calibration and validation t hrough a collaborative effort of academicians and practitioners in the area of H/WQ modeling (personal communication). PAGE 17 17 Restrepo and Schaake (2009) in their review on challenges and directions of research in hydrologic forecasting discuss current and futur e areas of research which include the coupling of surface and sub surface hydrologic models, combining physical and conceptual modeling approaches, and downscaling and upscaling issues while working with climate models. One aspect which is emphasized in th is article and other studies ( W agener and Kollat, 2007; Mirchi et al., 2009; Daniel et al., 2011 etc.) is the need to build confidence in models and their applications through reliability analysis. Uncertainty and sensitivity analyses of data and models are essential ingredients of modeling processes . They are also part of specific areas on the agenda of the ASABE task force committee on model calibration validation and a number of other recent efforts e.g. model uncertainty and sensitivity (Shin et al., 2013 ; Baroni and Tarantola, 2014 ), measurement uncertainty (Harmel et al . , 2006; Harmel and Smith, 2007), incorporation of both measurement and model uncertainty in model evaluation (Harmel et al., 2009; Ritter and Mu Carpena, 2013) etc . To summarize , recent research trends in H/WQ modeling suggest that considerable efforts are being done towards understanding and improving reliability of these models. Uncertainty and sensitivity analyses are two important tools included in the overall reliability eval uation framework (e.g. Refsgaard et al., 2007). Development of these techniques is still an ongoing process (e.g. Saltelli et al., 201 0 ; Kucherenko et al., 2012; Ge et al., 2014) . Also some recent works have suggested the need to standardize model evaluati on and reliability communication practices for H/WQ models ( e.g. Shin et al. 2013; Gan et al., 2014 ; Harmel et al., 2014 ) . Uncertainty and sensitivity analysis techniques and their applications in H/WQ modeling to achieve efficiency and understanding form the core of this research. PAGE 18 18 Sensitivity and Uncertainty Analysis: An Overview Sensitivity analysis (SA) can be defined as the study of apportioning the uncertainty in the model output to the uncertainty in the model input parameters while uncertainty analysis (UA) is the quantification of uncertainty in model output (Saltelli et al., 2004). Hence SA enables the modeler for (a) model corroboration i.e. validation or falsification of the model, (b) research prioritization i.e. detailed analysis or measur ement of influential parameters, (c) model simplification e.g. parameter fixing, (d) identification of regions in the parameter space that result in extreme outputs , (e) determination of optimal regions in the parameter space to be used in subsequent model calibration , and (f) if and which parameters interact with ea ch other (Saltelli et al., 2005; Saltelli et al., 2008). On the other hand , knowledge of predictive uncertainty of model output, as determined in UA, is of foremost concern when the model is bei ng used for any kind of decision making Carpena et al., 2006; Shirmohammadi et al., 2006) . Sensitivity Analysis Techniques Details of the various UA/SA methods can be found in Saltelli and Marivoet (1990), Melching (1995), Saltelli et al. (20 00), Helton and Davis (2002), Saltelli et al. (2004), Saltelli et al. (2005) , Manache and Melching (2008) , Gan et al. (2014) to name a few. SA techniques can be classified on three bases exploration of factor space , purpose of SA, and type of sensitivity measure (Figure 1 1). The most widely discussed c lassification is based on exploration of factor space . If only one factor is varied at a time (keeping all other factors fixed at base value ) in the vicinity of a base point then it is a local SA . On the ot her hand if all factors are varied simultaneously and the entire factor space is explored then it is known as a global SA. This characteristic of factor variation has important implications on the applicability of a SA technique and the interpretation of results . Based on the purpose for which SA is performed , SA techniques can be classified into (a) qualitative/ screening (mere identification of important and PAGE 19 19 unimportant factors) and (b) The third class ification is on the basis of type of sensitivity measure which range from simple scatterplots to complex variance decomposition based sensitivity indices . Figure 1 1 also shows SA techniques under each category. Details of some of the commonly used SA tech niques , their advantages, and disadvantages are presented below . One at a time One of the simplest and most commonly used sensitivity analysis methods is changing one factor at a time (OAT), while keeping all other factor s fixed at their base values ( Neari ng et al., 1990; Haan and Skaggs, 2003; White and Chaubey, 2005). Sensitivity is commonly measured by the ratio of the change in output to the change in input factor while all other factor s are held constant . Many argue that OAT i s a reasonable approach as any output change will unambiguously be due to the change in a single factor . However, this appro ach needs to be used cautiously for non linear models since their derivatives strongly depend on the point in the factor space at wh ich they are being evaluated (Wallach et al., 2006). Also, b approach does not take into account the interaction of factors (Saltelli et al., 2008) . Method of e lementary e ffects / Morris method The EE method can be considered an extension of the OAT approach. However, unlike OAT, the EE method is a global method since OAT experiments are repeated all over the factor space i n order to assess sensitivity. The first step in the EE method is to generate a factor sample which is done from t he k dimensional factor hyperspace divided into a p level grid . A p level grid implies that each factor can take only p distinct values {0, 1/( p 1), 2/( p r trajectories (or paths), each consisting of k +1 points are generated . Within any given trajectory two samples differ from each other in only one factor p /[2( p 1)]. The EE associated with the i th model factor is calculated as: PAGE 20 20 ( 1 1 ) wh ere y is the model output to the corresponding factor set , x i is the i th input factor . Since each trajectory provides one EE for each of the k factor s, r independent trajectories give r EE values for each factor which can be statistically analyzed. Morris (1991) proposed to use the mean ( i ) and standard deviation ( i ) calculated from the EE distribution associated with the i th factor as sensitivity measures, where i is an indicator of tot al or overall importance , and the standard deviation i indicates factor interactions and non linearity associated with the i th factor . The total computational cost associated with this method is r ( k +1) simulation s . Different articles report different numbers for r, ranging from 2 to 80 , though some studies (e.g. Herman et al., 2013 a ; Wainwright et al., 2014) have indicated that very few ( r <20) trajectories are sufficient for reliable segregation of important model parameters. Campolongo et al. ( 2007) proposed an improvement to this method in which the enhanced Morris measure is the mean ( * i ) of the distribution of the absolute values of the EE i . Campolongo et al. ( 2007 ) showed that * i is a better sensitivity measure than i for ranking factors in order of their overall importance . For some non additive, non monotonic models, EEs of a particular input may be positive or negative for different trajectories, which may cancel each other and thus not get reflected in the calculated i , ultimately not being identified as influential in sp ite of strongly affecting the model output. This enhanced Morris sensitivity measure has been widely accepted as the standard for screening exercises using the method of EE (e.g. Pujol, 2009; Campolongo et al., 2011). Campolongo et al. (2007) also proposed a modified sampling strategy to select the final r trajectories from a pool of about 500 to 1000 trajectories based on a distance maximization PAGE 21 21 criteri on . Pujol (200 9 ), Ruano et al. (2012) , Ge et al. (2014) , and Sohier et al. (2014) have also proposed diff erent types of parameter sampling strategies for the method of EEs. Correlation/regression methods A number of sensitivity analysis methods involve multiple random or quasi random model simulation s. Correlation or regression coefficients between factor s and model outputs are calculated and used as sensitivity measures. In the sensitivity analysis literature these methods are often referred to as correlation/regression methods or as sampling ba sed techniques. Typical sensitivity measures include the Pearso n product moment coefficient (PEAR), standardized regression coefficient (SRC), partial correlation coefficient (PCC), Spearman coefficient (SPEAR), standardized rank regression coefficient (SRRC) , and partial rank correlation coefficient (PRCC). Saltelli and Marivoet (1990), Helton and Davis (2002), and Manache and Melching (2008) have assessed the strengths and weaknesses of these methods where the lat t er have also summarized some additional correlation/regression based measures. PEAR is the usual linear correlation coefficient calculated on output variables Y i and input parameter X ij , while SRCs are calculated by fitting a multi linear regression model. The linear regression coefficients b j s obtained from this model are then multiplied by the ratio of var iance of the j th parameter to total model variance. The PCC between output Y and the input parameter X j is obtained from sequential regression models. Two models are constructed first: ( 1 2 ) ( 1 3 ) w here and are response variables of regression models. Next, new variables Y and are calculated. The correlation coefficient between these two new variables is the PCC between Y and X j . PAGE 22 22 One of the drawbacks of PEAR, SRC , and PCC is that they perform well only if the model is linear. Hence while applying these methods, it is recommended to calculate the coefficient of determination (R 2 ) for the regression model to test the degree o f linearity of the model. Saltelli et al. (2004) suggested R 2 > 0.7 as a limit of applicability for linear regression based sensitivity techniques. SPEAR, SRRC and PRCC are the rank based counterparts of PEAR, SRC and PCC, respectively, as they are calcula ted by the same procedure but on ranked inputs and outputs. Hence these can be used as sensitivity indices if the model is non linear but monotonic (continuous and smooth rising or falling response). SRC and PCC produce identical parameter rankings if the inputs are independent but may differ in case of correlated inputs. Iman and Helton (1988) proposed the use of PCC when parameters are correlated as they provide a measure of the strength of the linear relationship between output variables and parameters a fter a correction has been made for the linear effects of the other parameters. However, a rank based restrictive correlation structure proposed by Iman and Conover (1982) needs to be specified while sampling. The success of these methods depends on the ty pe of scheme used to sample the input factor s. Amongst several sampling techniques available, e.g., random sampling, importance sampling, and fixed samples, Latin hypercube sampling (LHS) has been found to produce stable results (Helton and Davis, 2002). I n LHS marginal distributions of factors are divided into intervals of equal probability. Then samples are drawn randomly from each interval for every parameter, which ensures the full coverage of the range of each variable (McKay et al., 1979). Regionalized s ensitivity a nalysis (RSA )/ Hornberger Spear Young m ethod The regionalized sensitivity analysis ( R SA) technique is a global method first presented by Spear and Hornberger (1980) and Hornberger and Spear (1981) for the purpose of Monte Carlo fi ltering . Based on model output values (generally some likelihood function) model runs PAGE 23 23 are pa rtitioned into behavioral (acceptable) and non behavior al (non acceptable) categories. Factor sets corresponding to behavioral and non behavioral categories are cla ssified accordingly as well. Distributions of a given factor from the two portions are compared to assess the statistical difference by the Kolmogorov Smirnov (KS) two sample test (Massey, 1951). If the difference is statistically significant then the mode l is considered to be sensitive to that factor . One advantage of this approach is that no assumption is made about the covariance of factor s, i.e., it is a non parametric sensitivity analysis method (Beven, 2001). However, R SA cannot identify the relative importance of a parameter within the behavioral class (Saltelli et al., 2000) making it s omewhat qualitative in nature. Multi o bjective regionalized sensitivity a nalysis (MO R SA) Multi objective generalized sensitivity analysis (MO R SA) was developed by Bast idas (1998) as an extension of R SA to make the method more suitable for models with a large number of factor s, which is typically the case of H/WQ models. In MO R SA the concept of a singular criterion for behavioral/non behavioral model runs is replaced by multi criteria filtering based on a Pareto optimal set. Essentially in MO R SA n behavioral/optimal factor sets are separated from m non behavioral/non optimal parameter sets based on N selection criteria, which are then tested factor . Variance b ased m ethods Variance based methods use analysis of variance (ANOVA) like techniques to decompo se the output variance into components. The main advantage of using these methods is that they work equally well for all types of models (e.g. linear, non linear, monotonic, non monotonic) (Saltelli et al., 2000). Two of the most commonly used variance bas ed global PAGE 24 24 A mplitude S ensitivity T est (FAST). function f ( X ) (i.e. model response Y ) into summations of increasing dimensionality as ( 1 4 ) evaluated via multidimensional integrals. The contribution of each term in eq. 1 4 to the total variance (unconditional variance) V can be easily obtained from ( 1 5 ) The equation can be standardized by dividing throug hout by V to give ( 1 6 ) where S i are the first order sensitivity indices, S ij the second order sensitivity indices, and so on, and k is the number of model factor s. The total effect index of the i th parameter is calculated by adding its first order index and its contribution to higher order indices. For example, if k = 3, then the total effect index for the 1 st factor is calculated as ST 1 = S 1 + S 12 + S 13 + S 123 S i = 1 indicates that the model is perfectly additive . In practice only the first order and total effect sensitivity indices are calculated. The calculations are done as follows: ( 1 7 ) PAGE 25 25 ( 1 8 ) where is the variance of expectation of Y over i th factor , is the va riance of expectation of Y over all but i th facto r and similarly is the expectation of variance of Y over all but i th factor . The computational cost for traditionally used algorithms to implement the is N (2 k +2), where N is the number of model runs per index. Saltelli (2002) extended the original the calculation of the first and the total effect sensitivity indices with N ( k +2) model runs. There are vario us variance decomposition schemes and two parameter sampling schemes; LP sequences and winding stairs that can be used for this method formula (Jansen, 1999) coupled with LP seq uences produces stable total sensitivity indices at lower N . Saltelli et al. (2005) recommended N in the range of 500 1000 to be sufficient for stable sensitivity indices. FAST and extended FAST (eFAST) : The Fourier amplitude sensitivity test (FAST) was first proposed by Cukier et al. (1973, 1978) and subsequently modified by a number of scientists. Originally FAST could only identify the first order indices. The main idea behind FAST is to convert multidimensional integrals (as dim ensional integrals in s by multidimensional Fourier transforms. A search curve is defined which explores the k dimensional unit hypercube parameter space by a set of parametric equations, (1 9 ) PAGE 26 26 where, x i is the i th parameter as a function of s , G i is the transformation function, s = [ scalar variable and i is a set of integer angular frequencies. Saltelli et al. (1999) extended this idea by defining a different transformation function, which is more efficient and allowed resampling by the introduction of a random phase lag. This enabled the calculation of total effects and the method is known as ex tended FAST or eFAST. Details of sensitivity indices (first and total order) for eFAST can be found in Saltelli et al. (1999). The cost of computation is N s k . N s is the number of model evaluations required per parameter which is estimated as N s ~ N r (2 max +1) where N r is the number of resamples (or phase lags), M is the interference factor and max is the maximum of parameter frequencies. Saltelli et al. (1999) also give a procedure to calculate frequencies assuming N r = 2 and M = 4. They recommend usin g N s ~ 500 1000 to get robust results. Recent advancements in SA techniques Some of the recent efforts in SA t echnique s have been focused on (a) improving the method of EE, especially the sampling strategies (e.g. Pujol, 2009 ; Campolongo et al., 2011; Ruan o et al., 2012; Ge et al., 2014; Sohier et al., 2014 ), (b) variance decomposition for the index (e.g. Wu et al., 2012 ; Wainwright et al., 2014 ) , (d) derivativ e based global sensitivity , and (e) SA techniques for dependent (correlated) parameters (e.g. Xu and Gertner, 2008; Kucherenko et al., 2012) . Uncertainty Analysis Techniques UA methods im plemented in H / WQ simulation studies range from output distribution statistics (e.g., mean, range, interquartile ranges), output histograms/distribution functions, first and higher order errors and mean estimates (e.g., Shirmohammadi et al., 2006) to quite sophisticated Markov Chain Monte Carlo (MCMC) methods. Stow et al. (2007) recommended PAGE 27 27 five UA methods for adaptive TMDL implementation, namely (1) RSA or Generalized Sensitivity Analysis, (2) Generalized Lik elihood Uncertainty Estimation (GLUE) , (3) Bayes ian Monte Carlo, (4) Importance Sampling, and (5) MCMC, all essentially Monte Carlo methods. Among these techniques RSA (Hornberger and Spear, 1981) is not really an UA method but forms the basis of the GLUE methodology (Beven and Binley, 1992), which is t he most widely used UA technique in H/WQ modeling . RSA classifies model runs into behavioral and non behavioral runs based on some fitting criteria. In GLUE, the concept of behavioral/non behavioral runs is replaced by weighting of factor sets based on likelihood measures of model runs. The likelihood measures used are generally goodness of fit criteria (e.g., root mean square error). Model runs can be classified into behavioral and non behavioral runs using some subjective likelihood crit eria and further rescaled so that the cumulative weight (of behavioral runs) is unity. However, there is considerable subjectivity involved at various stages in GLUE (e.g., prior distributions, likelihood function, etc.) which can hamper the results. MCMC employs a full Bayesian framework and uses unassisted algorithms to choose a sample that approaches the posterior density function. This makes MCMC the most efficient and fast method for UA. However, the efficiency of MCMC depends on the underlying algorit hms to locate the optimal factor combinations for the response surface. In H /WQ models response surfaces are quite complex; therefore, MCMC may need a large number of runs (Beven, 2001; Gallagher and Doherty, 2006). Ratto et al. (2001) proposed a G lobal S A GLUE approach for model calibration which respectively. This allows sensitivity analysis on raw model outputs as well as likelihood measures. The former helps in f actor fixing and research prioritization, while the later helps in PAGE 28 28 the parameter identification stage of model calibration. The GLUE part of this procedure estimates the predictive uncertainty. Since the same model runs are utilized for SA and UA, this is an efficient technique. Sensitivity and Uncertainty Analysis in H/WQ Modeling: Applications and Issues Details of a number of SA/UA applications (study, number of parameters, outputs etc.) in H/WQ modeling have been summarized in Tables 1 2 to 1 4. More elaborate summaries of SA and UA in H/WQ modeling can be found in Reusser et al., (2011), Yuan et al. (2014) , and Guzman et al. (2014). From these summaries it i s evident that most of the SA/UA applications were intended for model calibration and paramete r estimation as objective functions/ likelihood measures were considered as the model output s . Since the results of such applications are conditioned to the evaluation measure used, model sensitivity analysis results may not identify algorithms responsible for model behavior (Reusser et al., 2011). Also , in many applications it was assumed that parameters are the only source of uncertainty while model structure, assumptions, initial and boundary conditions , and input data can contribute to the overall uncertainty (Melching, 1995). The task of performing SA and hence calibr ation validation in the case of a large scale H/WQ model is t edious due to the complexity of the processes simulated by such models and t he issue of over parameterization (e.g. Cibin et al., 2010; Yang et al., 2012 ; Gan et al., 2014 ). The complexity of the problem further increases when the model is spatially distributed (Crosetto et al., 200 1; Tang et al., 2007 ; Gan et al., 2014 ). Spatiall y distributed models require significant quantities of spatially distributed input data, further adding to the parameter uncertainty (DeBarry and Quimpo, 1999). It is a substantial source of error that propagates through the model and affects the uncertain ty of results (Za j a c , 2010). Though some studies focused on PAGE 29 29 assessing the model uncertainty associated with the input data ( digital elevation models , rainfall, soil etc.) have been done (see Table 1 4 ), such sensitivity analys e s for input data are rare. F igure 1 2 and Table 1 5 indicate that in spite of a number of limitations of the OAT method (discussed earlier) it has been the most widely used technique with ~40% of overall applications. Some of the likely reasons behind this are ( a ) it is the least exp ensive method, ( b ) ignorance of its limitations by the modeling community, and ( c ) OAT has been implemented in built in SA tools of some watershed models (e.g. SWAT, Mike SHE) and parameter estimation packages such as PEST . The method of elementary effects i s the second most abundantly used methodology. Though designed for factor screening some of the studies have used it as quantitative technique (e.g. Srivastava et al., 2014) . are used less frequently in spite of their superiority to OAT and Morris approaches essentially due to the typical high dimensionality of H/WQ models requiring an impractically large number of model runs . Some studies have recommended the use of sequential SA i.e. formal factor screening through Morris analysis before performing quantitative SA using variance based global methods oz Carpena et al., 2007 etc.). However, many studies continue to selec t few parameters for variance based global SA based on previous model studies, conceptual understanding of the system, and/or expert knowledge. For the UA applications in H/WQ MCMC has been implemented in relatively few studies compared to GLUE (Table 1 3 and Guzman et al ., 2014 ) . Mishra (2011), in his work on evaluation of predictive uncertainty of HSPF in TMDL development using Monte Carlo based techniques, concluded that GLUE should be preferred over MCMC due to its flexibility and relative computational economy which b ecomes very important aspect in evaluation of watersh ed scale water quality models. The s uperiority of GLUE over MCMC or vice versa in H/WQ modeling is PAGE 30 30 still an unresolved issue. A c onsiderable amount of literature has been published in recent years by the proponents of these approaches (e.g. Vrugt et al., 2009; Beven et al., 2011 etc) . Another important aspect of a number of SA/UA techniques is the implicit assumption of orthogonality (i.e. absence of parameter correlation) of the input factor s e.g. OAT, FAST/eFAST etc. Among the several SA analysis techniques described previously only a few correlation/regression based techniques can take into account factor correlations e.g. PCC, PRCC. However, applicability of these methods is dependent on t he nature of the model (linearity and monotonicity) which is usually unknown. Considering the complexities of H/WQ models it would be tenuous to assume them to be linear and/or monotonic in order to apply correlation/regression based SA techniques (Gan et al., 2014) . In spite of recent development in variance decomposition based SA techniques for the correlated parameter case (Xu and Gertner, 2008; Kucherenko et al., 2012), this aspect of SA /UA applications in H/WQ modeling is not yet explored sufficiently . In fact none of the studies (based on Table 1 5) explicitly discussed and/or accounted for potential parameter dependencies. Also, as mentioned previously, high dimensionality is a typical issue associated with large scale H/WQ models. Reduction of overal l model uncertainty through qua litative segre gation of factors into important and unimportant categories (known as factor screening) before application of computationally intensive SA/UA techniques is recommended (e.g. Wang oz Carpena et al., 2007) . It was noted that the method of EE is the most widely used screening method and the recent EE technique development is focused on improving sampling strateg ies . Currently available sampling strategies are founded on either hyper distance based space coverage or do not use any selection criteria sample generation. Uniformity of the generated parameter sample is not handled explicitly. Besides, strategies proposed by PAGE 31 31 Campolongo et al. (2007) and Ruano et al. (2012) need impractically large sample generation time even when the model dimensionality i s moderately high. Hence, a new sampling strategy which will account for sample uniformity, space coverage, and will be computationally efficient is required. The second issue with the EE method is that ( a ) currently it is not possible to compare sensitivity indices across different outputs and ( b ) sensitivity indices are good only for parameter ranking . For the possible Campolongo et al. (2011) , it is necessary to improve EE sensitivity measures. In summary, evaluation of H/WQ models through uncertainty and sensitivity analyses is not n ew . However, the potential of these techniques as evaluation tools has been underutilized in practice . There are severa l areas of UA/SA implementations in H/WQ that need improvement which include (but not are limited to) (a) incorporation of factor correlation and improvement of evaluation tools for the same, (b) evaluation of H/WQ models for sensitivity to input data (cli matology, land use etc.), (c) improvement of available SA tools to handle high dimensionality of H/WQ models explicitly and efficiently , and (d) standardization of UA/SA practices in H/WQ. Objectives The overall goal of this work is to explore and improve sensitivity analysis techniques in non typical model evaluation settings for H/WQ models ( consideration of factor correlation and high dimensionality) in order to further increase the value of these models as reliable decision making tools. The Agricultural Reference Index for Drought (ARID) a drought index model for southeastern United States and the Watershed Assessment Model (WAM) a spatially distributed watershed scale model will be used as test cases for sensitivity and uncertainty analy sis applications. T he specific objectives of this work are PAGE 32 32 T o perform global sensitivity and uncertainty analysis of a simple drought index model with correlated parameters using two different sensitivity analysis techniques . T o develop and evaluate a new multi criteria trajectory based parameter sampling strategy for the screening type sensitivity analysis method of Elementary Effects , an important step in SA of high dimensional models . T o implement the newly developed parameter sampling strategy as a par t of parametric sensitivity and predictive uncertainty model evaluation framework for a high dimensional watershed scale model . A study based on the first of the three objectives is presented in Chapter 2. The development of new smapling strategy and its c omparison with other benchmark strategies is explained in Chapter 3. Chapter 4 presents the evaluation of Watershed Assessment Model, study corresponding to objective three. The last chapter i.e. Chapter 5 summarizes the major findings, limitations, and fu ture research needs of this work. PAGE 33 33 Table 1 1. Summary of widely used watershed scale models 1 Model Components Spatial Scale and Watershed Representation Temporal Scale Availability and Application History Suitability AnnAGNPS (Bingner and Theurer, 200 1 ) Hydrology, sediment, nutrients and pesticide transport Distributed, Homogenous land areas, reaches and impoundments Continuous, Daily or sub daily time steps Public. Supported by USDA. Predominantly used in U.S. Agricultural watersheds for the eval uation of management practices ANSWERS 2000 (Bouraoui and Dillaha, 200 2 ) Runoff, infiltration, river routing, chemical/ nutrient and sediment transport Distributed. Square grids/ cells inside which soil/ land use are homogenous Continuous, Variable time steps Public. Predominantly used in U.S. Medium size agricultural watersheds, evaluation of BMPs, capable of handling transformation and interaction between four nitrogen pools HSPF (Bicknell et al., 2001) Runoff, water quality constituents, pervious and impervious surfaces, stream channels Semi distributed. Pervious/ impervious land area, stream channels and mixed reservoirs, 1 D simulation Continuous. Variable but constant time steps (hourly) Public Widel y used in U.S. and other parts of the world. Supported by USEPA and USGS Suited for both agricultural and urban watersheds, diverse water quality and sediment transport MIKE SHE ( Refsgaard and Storm, 2012 ) Interception, overland and channel flow, unsaturated and saturated zone flow, aquifer river exchanges, advection/dispersion, geochemical processes, soil erosion Distributed. Square grids for overland flow, 1 D channels and unsaturated zone and 3 D saturated zone Continuous. Long term and event ba sed. Variable time steps Private. Has been applied in various conditions all over the world Wide range of spatial and temporal scales, modular design, parameter estimation and water budget analysis SWAT (Neitsch et al., 2002) Hydrology, weather, sediment , soil temp., crop growth, nutrients, pesticides, channel and reservoir routing Semi distributed. Sub watersheds based on climate, HRU based on soil, cover and management, ponds, groundwater and main channels Continuous. Daily time steps Public. Supported by USDA and USEPA. Extensively applied all over the world Best suited for agricultural watersheds, excellent for TMDL calculations and BMP evaluations. WAM (Bottcher et al., 2012) Runoff, water quality constituents, pervious and impervious surfaces, stream channels Distributed. Square grids/ cells inside which soil/ land use are homogenous Continuous. Daily time steps Public. One of the US EPA recommended model. Extensive applications in the state of Florida Can efficiently handle watershed with diver se land uses a nd hydrologic conditions. TMDL and BMP evaluation 1 Based on Borah and Bera (2003), Migiaccio and Srivastava (2007) , and Daniels et al., (2011) PAGE 34 34 Table 1 2. Summary of selected sensitivity analysis applications in H/WQ modeling studies Study Model Parameters Method 1 Runs Output 2 Deflandre et al. (2006) QUESTOR 5 12 parameters eFAST 500 NSE Cloke et al. (2008) ESTEL 2D 9 parameters M OR SA <1280 Fuzzy membership, sum squared errors van Werkhoven et al. (2009) SAC SMA 14 parameters for upper and lower zones, partition and percolations 7.5e6 4 objective functions van Werkhoven et al. (2008 ) SAC SMA 78 cells x 14 parameters 4e6 RMSE Wagener et al (2009) SAC SMA 78 cells x 14 parameters Not reported RMSE, TRMSE, ROCE Tang et al. (2007) SNOW 17 SAC SMA 5 parameters for snow and 13 for SAC SMA Sobol, RSA, ANOVA, LSA <10000 RMSE of discharge, RMSE of Box Cox transformed discharge Pappenberger et al. (2008) HEC RAS 6 parameters Sobol, Morris, RSA and Regression Not reported Mean abs error, NSE Pappenberger et al. (2006) HEC RAS 7 parameters Regional split and regression trees, correlation, RSA 3000 Inundation performance measure Demaria et al. (2007) VIC 10 parameters RSA, scatterplots 60000 RMSE, RMSE Box Cox, ARE of discharge McIntyre et al. (2003) WaterRAT 26 parameters RSA 10000 Water quality concentrations Sieber and Uhlenbrook (2005) TAC 20 parameters RSA, regression 400 Discharge, NS E of discharge, Christiaens (2002) MIKE SHE Soil hydraulic parameters Regression 25 multiple Muleta and Nicklow (2005) SWAT 35 parameters Regression 300 RMSE Van Griensven et al. (2006) SWAT 41 parameters Modified Morris, hybrid of OAT and LHS Not reported Total sum of squared errors of flow, sediment and nutrient Foglia et al. (2009) TOPKAPI 35 parameters Local SA 71 discharge Benke et al. (2008) 2C Model 5 parameters Stepwise local SA 30000 Annual discharge PAGE 35 35 Table 1 2. Continued Study Model Parameters Method 1 Runs Output 2 Cullamann et al. (2006) WaSiM ETH 6 parameters Loacal SA 13 Peak discharge Nossent et al. (2011) SWAT 26 parameters bootstrapping 336000 NSE Arabi et al. (2007) SWAT 41 parameters MO R SA, TSDE 5000 NSE Bahremand and De Smedt (2008) WetSpa 11 parameters PEST Not reported NSE Francos et al. (2003) SWAT 22 parameters Morris, FAST Not reported Multiple Smith and Wheater (2004) P export and discharge 8 parameters RSA 30000 RMSE of total discharge and TP concentrations Vezzaro and Mikkelsen (2012) Micro pollutants in streamflow 10 parameters 10000 Multiple Yang et al. (2012) WetSpa 16 parameters Morris Method and State Dependent Parameter Method <850 NSE, MAED, MAEC Reusser et al. (2011) TOPMODEL and WaSiM ETH 9 parameters for TOPMODEL and 11 parameters for WaSiM ETH FAST, eFAST and <5632 Discharge time series Zheng and Keller (2006) WARMF 121 RSA 10000 Multiple Wang et al. (2006) APEX 15 Morris and FAST Not reported Shen et al. (2008) SWAT 21 Morris <100 Multiple Carpena et al. (2007) VFSMOD 27 Morris and eFAST Not reported Multiple Herman et al. (2013) SAC SMA 1093 Morris and Sobol 21860 6.5e6 RMSE related to streamflow outputs 1 MORSA Multi objective regionalized sensitivity analysis, FAST Fourier amplitude sensitivity test, RSA Regionalized sensitivity analysis, ANOVA Analysis of variance, PEST parameter estimation and sensitivity test, OAT one at a time, LSA local sensitivity analysis, LHS L atin hypercube sampling 2 NSE Nash Sutcliff e efficiency, RMSE r oot mean square error, TRMSE Transformed RMSE, ROCE runoff coefficient error , ARE absolute relative error, MAED mean absolute error for duration curves, MAEC mean absolute error for accumulation curve PAGE 36 36 Table 1 3. Sum mary of parametric uncertainty analysis applications in selected H/WQ modeling studies Study Model Parameters Method 1 Runs Evaluated Output 2 O sidelo et al. 2003 Sediment nutrient dynamic model 5 management and 14 sediment and nutrient related parameters MC technique: RSA UCPR TSDE Not reported Multiple criteria Gallagher and Doherty (2007) HSPF 7 MCMC, FOEA >250 -Beven and Binley (1992) IHDM 4 GLUE 500 (x 5 storm events) NSE, sum of squares of the residuals, scaled maximum absolute error Migliaccio and Chaubey (2008) SWAT 13 parameters, 12 HRU configurations Single factor between subjects analysis 500(x 12 HRU combo) % relative error, coefficient of variability Muleta and Nicklow (2005) SWAT 35 GLUE 5000 NSEs for discharge and sediment concentration and their combinations Smith and Wheater (2004) P export and discharge 8 parameters GLUE 30000 RMSE of TP concentrations Vezzaro and Mikkelsen (2012) Micro pollutants in streamflo w 5 most sensitive parameters from SA GLUE 10000 Multiple Wu and Liu (2012) SWAT R 8 MCMC 1000 Multiple criteria on streamflow Yu et al. (200 1 ) Distributed rainfall runoff model 3 (5 storm events) MCS, LHS, RPEM, HPEM <1000 Different for different methods Zhang et al. (2009) SWAT 11 Bayesian Model Averaging 1000 Wang et al. (2006) DRAINMOD 8 GLUE 1992 MSE and its variants for flow Shen et al. (2008) SWAT 21 MC, FOEA <100 Multiple Vrugt et al. (2008) HYMOD 5 MCMC -Multiple Blasone et al. (2008) HYMOD, NAM and SAC SMA 5 for HYMOD, 10 for NAM and 14 for SAC SMA Combination of MCMC and GLUE -Exponential likelihood function Jin et al. (2010) WASMOD 3 GLUE and Bayesian method 50000 Multiple based on flow Carpena et al. (2007) VFSMOD 27 CDFs and other statis t i cs of outputs <8000 Multiple 1 UCPR Uniform covering by probability rejection, TSDE tree structured density estimation, MCMC Markov Chain Monte Carlo, FOEA first order error analysis, GLUE generalized likelihood uncertainty estimation, MCS Monte Carlo simulations, LHS La tin hypercube sampling, RPEM point estimation method, HPEM 2 NSE Nash Sutcliff e Efficiency, RMSE root mean square error PAGE 37 37 Table 1 4. Summary of input data uncertainty analysis applications in selected H/WQ modeling studies Study Model Analysis Details Chaubey et al. (2005) SWAT Effect of DEM resolution on output uncertainty 7 DEM resolution scenarios Chaubey et al. (1999) AGNPS Effect of spatial rainfall variability on parameter uncertainty 17 rainfall distribution scenarios Cho et al. (2009) SWAT Effect of spatial distribution of rainfall on temporal and spatial uncertainty (1) Effect of rainfall input method in conjunction with subwatershed delineation (9 scenarios) (2) Effect of raingauge density on hydrology and water quality (10 scenarios) Cotter et al. (2003) SWAT Effect of spatial resolution of DEM, land use and soil on water quality model output uncertainty 7 resolution scenarios Muleta et al. (2007) SWAT Sensitivity of SWAT to spatial scale Calibrated parameter values for six different discretization scenarios were compared Sudheer et al. (2007) SWAT Impact of time scale of objective function (e.g. daily, annual etc.) on model performance Model perfo rmance measures were compared by for daily and monthly time scales Tang et al. (2007) SAC SMA Parameter sensitivity analysis at different time scales Four different SA methods used to evaluate effect of (3 different) time scales on parameter sensitivity van Werkhoven et al. (2008) HL DHMS Effect of rainfall characteristics and initial moisture conditions on spatial parameter sensitivity Total 6 combinations of rainfall and initial soil moisture characteristics were used spatial SA. For each case spatial SA indices were calculated. Miller et al. (2007) KINROS2 Effect of land use misclassification on output (runoff) uncertainty 100 realizations based on land use misclassification error were generated Patil et al. (2011) HSPF Effect of temporal and spatial data resolution on output uncertainty (1) Climate data (temperature and rainfall) at 5 different temporal resolutions were used to assess the impact on model output (2) 5 different resolutions of land use data were used to assess the impacts of spatial resolution Wu et al. (2005) Empirical Soil Loss Model Effect of DEM resolution on soil loss 7 different DEM resolutions were utilized Wu et al. (2007) TOPMODEL Effect o f DEM resolution and elevation data uncertainty 7 DEM resolution scenarios and thousands of realizations based on elevation errors are analyzed PAGE 38 38 Table 1 5. Summary of applications of sensitivity analysis techniques in H/WQ modeling pre and post 2000 based on Table 1 2 and Yuan et al. (2014) SA Technique 1 Number of Studies pre 2000 post 2000 Total OAT/LSA 16 30 46 EE/Morris 18 18 Correlation/Regression 1 8 9 RSA/MO R SA 10 10 FAST/eFAST 11 11 Sobol' 14 14 Other s 2 5 7 1 OAT One factor/parameter at a time, LSA local sensitivity analysis, EE elementary effect, RSA regionalized sensitivity analysis, MORSA multi objective regionalized sensitivity analysis, FAST Fourier amplitude sensitivity test, eFAST exten d ed Fourier am plitude sen sitivity test PAGE 39 39 Figure 1 1. Classification of sensitivity analysis techniques PAGE 40 40 Figure 1 2. Relative usage of various sensitivity analysis techniques applied in H/WQ studies. Results are based on 115 studies (Table 1 2 and Yuan et al. 2014) OAT/LSA 40% EE/Morris 16% Correlation/R egression Based 8% RSA/MOGSA 9% FAST/eFAST 9% Sobol' 12% Other 6% PAGE 41 41 CHAPTER 2 PARAMETER VARIABILITY AND DROUGHT MODELS A STUDY USING THE AGRICULTURAL INDEX FOR DROUGHT (ARID) A gricultural systems, especially rain fed and in arid regions, are prone to drought which ultimately affect s crop yield. Drought indices such as the Crop Moisture Index (CMI) (Palmer, 1968), the Soil Moisture Deficit Index (SMDI) (Narasimhan and Srinivasan, 2005), the Standardized Precipitation Index (SPI) (McKee et al., 1993), and the Water Requirement Satisfac tion Index (WRSI) (Frere and Popov, 1986) have been developed to quantify the severity of water shortage . While most of these indices consider the soil plant atmosphere continuum in their formulation, their applicability suffers as : (1) some are crop speci fic indices ( e.g. SMDI and WRSI), requiring a large number of input parameters ; (2) some do not consider soil and/or plant characteristics , (e.g. SPI and the Crop Water Stress Index (CWSI) (Idso et al., 1981)); and (3) some employ large time steps (weekly, monthly , etc.) that are not optimal for estimating crop water stress, such as the SPI and the Palmer Drought Severity Index (PDSI) (Palmer, 1965). A comprehensive summary of drought indices and their characteristics can be found in Mishra and Singh (2010) and Woli et al. (201 2 ) . In the latter, they have further developed a new generic daily water balance drought index called ARID. It has been found to simulate crop water balance dynamics fairly well for a number of crops, soils and weather conditions in th e southeastern USA (Woli et al., 2012; Woli et al., 2013). Further it has been implemented as a tool to monitor drought at various locations in the southeastern USA (Fraisse et al., 2011). Water stress forms the central part of a number of crop models suc h as the Decision Support System for Agrotechnology Transfer model (DSSAT; Jones et al., 2003), Aqua Crop (Raes et al., 2009), and the FAO56 water balance approach procedure (Allen et al., 1998); hence their successful application requires careful determin ation of water stress. Mathematical model s of any kind are abstraction s of reality due to incomplete understanding of governing physical PAGE 42 42 processes (physically based models), assumptions and approximations , data used for model development (empirical models) , model structure, measurement and instrumental errors , etc. ; all of which introduce uncertainties of various types ( e.g. Beven, 1989 ; Shirmohammadi et al., 2006 ). Sensitivity analysis (SA) can be used to quantify the strength of the input output relations hips of mathematical models while uncertainty analysis (UA) can be used to evaluate the propagation of uncertainties in input s onto model outputs (Saltelli et al. , 2004). UA/SA is a fundamental tool for understanding model behavior, the potential for model simplification, the identification of important input parameters and their probability distributions (e.g. Monte Carlo filtering) to aid model calibration process and reducing predictive uncertainty, and model development as a whole ( Tarantola and Saltell i, 2003 ; Saltelli et al., 2005 ) and is thus essential for the successful application of any mode l (Confalonieri et al., 2010). Neglecting to conduct UA/SA could result in undermining the science and value of the model (Beven, 2006). Methods of SA can be broadly classified into two groups , local and global. In local SA (LSA) methods , the local derivative of the desired output variable is numerically evaluated around a certain value of one input parameter , holding other input parameters constant at their me an values ; while in global SA (GSA) methods all parameters are varied over the entire parameter space (Saltelli et al., 2008) and are thus preferred over LSA methods given their ability to capture both the direct effects and the interaction effects. GSA me thods can be further classified into three types: (1) s creening methods , such as the elementary effect method of Morris ( Morris, 1991; Campolongo et al., 200 7); (2) correlation and regression indices such as Pearson product moment coefficients and Spearman rank correlation coefficients, and their partial and ranked counterparts; and (3) v ariance based methods (ANOVA like techniques), such as the Fourier Amplitude Sensitivity Test (FAST) (Cukier et al., 1973; Cukier et al., 1978) , PAGE 43 43 extended FAST (Saltelli et al., 1999) , and the method of Sobol GSA methods, especially variance based methods, are computationally expensive and as a result their applicability becomes an issue when the model has many parameters (Saltelli et al., 2005). Variance bas ed methods proportion the total output variance to individual parameters and their interactions. Unlike correlation and regression based methods, variance based methods do not suffer in their applicability due to model non linearity or non monotonicity (Sa ltelli et al., 2000; Dimov and Geor gieva , 2010) . One of the issues with variance based methods is that they can only be used when the input parameters are orthogonal (i.e. non correlated parameters). However, some of the correlation/regression based metho ds are capable of dealing with non orthogonal parameters. There are very few studies in the environmental and ecological sciences which deal with the parameter correlation issue. Some of the approaches proposed to tackle this situation include partial corr elation coefficients (PCCs) (Iman et al., 1985), Gram Schmidt orthogonalization (Bedford, 1998), correlation ratios (Saltelli et al., 2004), and decomposition of parameter variation into correlated and independent components (Xu and Gertner, 2008). The method of PCCs has been discussed in a number of works (e.g. Iman and Helton, 1988; Whiting et al., 1993; Iman et al., 2002) and has been used with some success. Prats and PicÃ³ (2010), Salemi et al. (2011), Padilla et al. (2011) , Woli et al. (2013), and o thers have performed sensitivity analyses of drought indices or water balance models. However, these studies either employed LSA techniques or only considered a single study site or only considered inter soil class variability or combinations of these. Thi s article aims at UA and SA of ARID considering inter soil and intra soil class parameter variability, plant characteristics, and location impacts to identify the most sensitive and influential model inputs along with nature of PAGE 44 44 ARID by applying two types o f global SA and UA techniques. Analyses were performed for five locations in the southeastern USA for four different soil types to assess the effect of location/climatic conditions. As soil hydraulic properties are typically correlated, correlation based s applied by accounting for the expected correlation of soil hydraulic properties in order to account for this shortcoming in variance based methods. Results of the two correla tion based SA M aterials and M ethods Model The ARID is based on a simple soil water bal ance that is similar to the FAO 56 approach (Allen et al., 1998) used for irrigation scheduling at the farm level. It u ses a reference crop approach making it independent of dynamic and crop specific coefficients. It considers the soil plant atmosphere continuum and operates at a daily time step. Due to these features ARID is ex model (Woli et al., 2012; Woli et al., 2013). Plant water stress can be expressed as the ratio of water deficit to water need (Thornthwaite, 1948), based on which the ARID is calculated as: ( 2 1) where ET = actual evapotranspiration (mm) and ET o = reference evapotranspiration (mm). When there is adequate moisture in the soil so that ET is equal to ET o , the value of ARID is zero (no stress), when soil moisture reaches the wilting point the value of ARID is unity (maximum stress). In maintaining ARID as a generic index, it does not employ crop coefficients to modify ET o , though it should be acknowledged that evapotranspiration from turfgrass surfaces in the PAGE 45 45 southeastern USA can differ markedly from refere nce ET o estimated for a hypothetical reference grass (Carrow, 1995; Romero and Dukes, 2009). The water balance components include precipitation, surface runoff, deep drainage, and ev apotranspiration (Figure 2 1). Capillary rise in the presence of a shallow water table is not included. The soil water content on the i th day is calculated from the soil water content on the (2 2) where, W i = soil water (mm) on the i t h day, W i 1 = soil water (mm) on the i 1 th day, P i = precipitation (mm) on the i th day, D i = deep drainage (mm) on the i th day, and R i = surface runoff (mm) on the i th day. Reference e vapotranspiration (ETo) Reference evapotranspiration is calculated using the FAO 56 Penman Monteith method (Allen et al., 1998). It assumes a hypothetical reference grass surface with crop height of 0.12 m, a fixed surface resistance of 70 s m 1 and an albedo of 0.23: (2 3) where R n = the net radiation at the crop surface (MJ m 2 day 1 ), G = the soil heat flux density (MJ m 2 day 1 ) (assumed to be negligible for daily application (Allen et al. 1998)), T = mean daily air temperature at 2 m height (Â°C), u 2 = the wind speed at 2 m heigh t (m s 1 ), e s = saturation vapor pressure (kPa), e a = t he actual vapor pressure (kPa), = the slope of the vapor pressure curve (kPa Â°C 1 ), and = the psychrometric constant (kPa Â°C 1 ). PAGE 46 46 Evapotranspiration (ET) ARID takes a macroscopic water content based approach to model actual evapotranspiration. It is based on the assumption that perennial crops transpire all water absorbed by roots from the soil (Wu et al., 1999) and roots are evenly distributed. It follows the water uptake model developed by Meinke e t al. (1993) using the exponential decay relationship proposed by Passioura (1983) (modified from Woli et al., 2012) : (2 4) where: = the volumetric soil moisture content (m 3 m 3 ), wp = the volumetric soil moisture content at the wilting point (m 3 m 3 ), MUF = the maximum uptake factor (day 1 ), and Z = the root zone depth (mm). Surface r unoff (R) The NRCS curve number method was adopted to calculate surface runoff (USDA, 1986) due to its simplicity and the availability of the required parameters. According to this method runoff is generated only when precipitation exceeds the initial abstraction: (2 5) where: CN = the runoff curve number and is dependent on soil type and land cover, I = the initial abstraction (mm), and S = the potential maximum retention (mm). PAGE 47 47 Deep d rainage (D) Soil can hold water only up to its field capacity and excess water is drained to lowe r soil layers below the root zone. In ARID deep drainage is calculated as (modified from Woli et al., 2012): ( 2 6 ) where: fc = the volumetric soil moisture content at field capacity (m 3 m 3 ) and DC = the drainage coefficient (m 3 m 3 day 1 ). Application/ Case Study Details Drought conditions at five locations in the southeastern USA, which included Smithfield, North Carolina (35.52 Â°N, 78.35 Â°W), Gainesville, Georgia (34.30 Â°N, 83.86 Â°W), Tifton, Ge orgia (31.45 Â°N, 83.48 Â°W), Saint Leo, Florida (28.34 Â°N, 82.26 Â°W), and Perrine, Florida (25.58 Â°N, 80.44 Â°W) (Figure 2 2), were simulated for a 24 year period using ARID. These locations were chosen due to their proximity to high quality weather stat ions contained within the United States Historical Climatology Network (USHCN) (Menne et al., 2012). These locations generally span the range of climate in the southeastern USA, from shorter growing seasons in the north to locations with distinct wet and d ry seasons in the case of the two Florida locations (Figure 2 3). Additionally four soil classes were considered. The effects of inter and intra soil class variability, location (climate) and variability in other model parameters on essed using two types of GSA/UA methods. GSA/UA Procedure A step by step procedure adopted for performing global sensitivity analysis include d: (1) determination of probability distribution functions (PDFs) of input parameters , (2) generation of PAGE 48 48 input samp les based on PDFs and the choice of SA method, (3) model simulations to calculate desired outputs, and (4) statistical analysis to obtain sensitivity ind ices and parameter rankings. SimLab v2.2 (Saltelli et al., 2004) developed at the Joint Research Centre of the European Commission, Ispra, Italy was used for pre processing (sample generation) and post processing (statistical treatment to determine sensitivity indices), i.e. step 2 and 4 respectively. The sample file generated from SimLab was used as the in put to the model which was simulated outside SimLab environment. O utput obtained from the model runs was stored in the desired SimLab post processing format for further analysis. Simulation S et U p Input data and p arameter PDFs Meteorological d ata. Daily precipitation, incoming solar radiation, maximum temperature, and minimum temperature for the period 1980 2003 for the five locations were obtained from the Daymet data center ( http://daymet.ornl.gov/ ) (Thornton et al., 1997; Thornton an d Running, 1999) (Figure 2 2). Daymet is a model developed by the Numerical Terradynamic Simulation Group at the University of Montana that produces serially complete estimates of daily weather variables at a 1 km resolution based on available ground based measurements. The Daymet dataset is constructed using an interpolation of observed precipitation and temperature along with their relationship to elevation us ing a digital elevation model. Solar radiation is estimated using diurna l temperature range and Sun Earth geometry following the methods of Bristow and Campbell (1984). Daily dew point temperatures were assumed to be equivalent to minimum temperatures, a valid assumption for reference conditions in humid regions (Allen et al. , 1998). Meteorological data were treated as deterministic inputs. P arameter d istributions. The two plant characteristics that influence ARID are Z and MUF . While the reference crop in ARID is a hypothetical grass, it resembles properties of perennial PAGE 49 49 turf grass with complete ground cover and uniformly distributed roots (Woli et al., 2012). Hence, the properties of turf grass were assumed to be applicable. Grasses have variable rooting depths which depend on their type, soil layering, and type of water mana gement. Grasses growing in layered soils and under irrigation tend to have shallower root depths while those growing in single layered soils and are rain fed tend to have deeper root depths. Cool season turf grasses have root depths ranging from 150 mm to 450 mm while warm season turf grass roots typically grow between 450 mm to 1500 mm ( Wu, 1985; Harivandi et al., 2009 ). Allen et al. (1998) recommended a general range of 500 1000 mm for the maximum root zone depth; further suggesting the use of smaller val ues of Z for irrigation scheduling and larger values for modeling soil water stress or rainfed conditions. Woli et al. (2012) used 400 mm as the default value for the root zone depth in ARID. Considering the wide range of observed and recommended values of root zone depths, following Woli et al. (2013) a uniform distribution was assumed for Z but with a wider range of 250 800 mm. While the MUF has been categorized here as a plant characteristic, it is also dependent on soil type (Dardanelli et al., 1997). Dardanelli et al. (2004) derived a generic value of 0.096 day 1 for MUF which was applicable to a number of crops and soil types. This value for MUF was used as the default value by Woli et al. (2012). However other studies (Meinke et al., 1993; Robertso n et al., 1993 a,b ; Dardanelli et al., 1997 ) have shown the value of MUF to vary from as low as 0.05 day 1 to as high as 0.12 day 1 . Based on data from these studies, MUF was fitted to a normal distribution with a mean 0.094 day 1 and standard deviation 0.0061 day 1 , for all soil types (Table 2 1). The soil hydraulic properties required by ARID are mainly governed by soil particle size distribution. For simplicity of analysis, four textural soil groups: Sandy Loam (SL), Silty Loam PAGE 50 50 (SiL), Sandy Clay Loam (SCL), and Silty Clay (SiC), based on the United States Department of Agricultural (USDA) textural soil classification triangle, were selected ensuring the representation of all four hydrologic soil groups defined in the Urban Hydro logy for Small Watersheds manual (USDA, 1986). The required soil hydraulic/hydrologic properties are wp , fc , CN , and DC (Table 2 1). The runoff curve number ( CN ) is a function of hydrologic soil group and type of cover (e.g. impervious, vegetated, etc). SL, SiL, SCL and SiC belong to hydrologic soil groups A, B, C, and D respectively. In the absence of well defined, continuous relationships for CN based on soil texture properties, ranges of CN were derived from values in TR 55 (USDA, 1986), and uniform probability distributions were assigned (Table 2 1). For DC , ARID uses a generic value of 0.55 assuming medium textured soils (Woli et al., 2012). For this study PDFs of DC were fitted based on the soil characteristic database of the DSSAT model (Jones et al., 2003) with modifications based on Suleiman and Ritchie (2004). A triangular distribution was fitted for SL while for the other soil types uniform distributions were assigned (Table 2 1). Texture based pedotranfer functions (PTFs) have been developed f or different regions by many researchers (e.g. Saxton et al., 1986; Rawls and Brakensiek, 1989; WÃ¶sten et al., 1999; Saxton and Rawls, 2006) to determine soil hydraulic properties such as fc and wp . These PTFs and their geographic origin have been noted to substantially influence model outputs (Loosvelt et al., 2011). Upon systematic comparison of various methods for estimating soil water retention parameters for crop models, Gijsman et al. (2003) concluded that the method of Saxton et al. (1986), which r equires the fraction of sand and clay, performs the best. Saxton and Rawls (2006) enhanced the PTFs of Saxton et al. (1986) based on updates of the USDA/NRCS National Soil Characterization database (Soil Survey Staff, 2004). PAGE 51 51 Values of wp and fc were calculated using the Saxton and Rawls (2006) PTFs, for a large number of points (at least 2000) uniformly sampled within the %sand %clay space for each of the four soil classes. Probability distribution functions were then fitted to the distributions of wp and fc for each soil class (Table 2 1). The wp and fc showed a strong positive correlation (correlation coefficient ~ 0.75). While the correlation between wp and fc was expected, it remains unclear to what extent the correlation values obtaine d can be attributed to actual soil data (physical effect) on which the PTFs were based or to the PTFs themselves (a mathematical artifact). The degree to which the correlation between wp and fc was caused by the latter was not accounted for in this analy sis. Due to the non orthogonality of wp and fc directly. Instead of using wp and fc as separate parameters they were combined into a single SHP ). A look up table with 100 equi probable points in soil space and corresponding wp and fc values was generated for each soil type. Each of the above points was randomly assigned a unique identification number (from 1 to 100). Hence SHP essentially has a discrete uniform probability distribution (Table 2 1). During sample generation SHP was sampled along with Z , CN , MUF , and DC . Values of wp and fc corresponding to SHP GSA/UA Methods The next step in the statistical framework is parameter sample generation for model executions. For this study we employed three sensitivity analysis methods namely, PCCs, ion based methods and are capable of analyzing model sensitivity in the presence of parameter correlation. While the method of PCCs suits well for linear and monotonic models, the method PAGE 52 52 of PRCCs is suitable for non linear but monotonic models. Since the n ature of the ARID model (linearity or monotonicity) was unknown both method types were used. For the correlation based analyses using PCCs and PRCCs, parameter samples were obtained using the latin hypercube sampling technique (LHS) (McKay et al., 1979) . LHS has been found to produce robust results with as few as 100 200 samples (Iman and Helton, 1988) compared to other methods such as random sampling, fixed sampling etc. (Helton and Davis, 2002) which require a much higher number of samples. However, it i s recommended to run the model as many times as possible for all of the above methods. Considering the very short computational time required for a single run of ARID, 10,000 samples were generated. To account for the correlation between wp and fc , restr icted pairing based on rank correlation structure (Iman a nd Conover, 1982) was defined. Though the correlation structure is based on ranks, the actual sample correlation approximately equals the desired values when the distributions are not highly skewed ( S i ) and total effect ( ST i ) sensitivity indices (which include first and higher order effects, where the higher order effects measure interactions) were calculated using k = 5 (number of model parameter s) and a sample size of N (2 k +2) = 6144 for each soil type and location combination, where N is the number of runs per model parameter. It is recommended to set N greater than 500 (Saltelli et al., 2004; Saltelli et al., 2008) and hence N = 512 was used for this study. Model Output, Uncertainty and Sensitivity Indicators Using the ARID index, Fraisse et al. (2011) classified drought threat into three zones depending on the value of ARID days in the high threat range (HT) was considered as the output variable of interest in this work. There are various ways to map the uncertainty between model parameters and outputs. Using only means and variances is not a very useful way of summarizing uncertainty a s there is PAGE 53 53 always some loss of information (Saltelli et al., 2000). Ibrekk and Morgan (1987) recommend displaying cumulative distribution functions (CDFs), density functions, and the mean value in the same figure while mapping uncertainty. In this study hi stograms of HT along with corresponding CDFs were plotted apart from summarizing output statistics. Scatterplots are very useful in identifying the nature of relationships between inputs and outputs (e.g. linear, monotonic, threshold, etc.) (Helton et al., 1993; Hamby, 1994) and can be regarded as qualitative and visual estimate of sensitivity. Scatterplots of HT against all parameters (for correlation analysis) were generated for all soil type and location combinations. Sensitivity indices obtained in the three methods are reported along with corresponding parameter rankings. Also coefficient of model determination i.e. R 2 adj values for raw as well as method, pa rameter rankings are based on total effects. Results Histograms of HT and other summary statistics, Figure 2 4 and Table 2 2, indicate the variability of the output at all five locations a nd four soil groups. For a given soil type, HT increased from north to south (from Smithfield to Perrine). Distributions were positively skewed for all soils with SL and SCL showing greater skewness compared to SiL and SiC. Central tendencies (mean values) of HT for SL and SCL showed similar values (though not identical) with values ranging from 0.11 to 0.23. However, the variability in SL, as seen from 95% confidence interval (CI) limits (Table 2 2), was greater than SCL and was the greatest of all soil ty pes. Silty loam had the lowest distribution means (0.06 to 0.14), but CI widths were similar to those of SiC for the corresponding locations. PAGE 54 54 Scatterplots (based on model runs for correlation analysis) Figure 2 5 shows the scatterplots of HT against all s ix parameters (partial correlation analysis) for the four soil classes under consideration. HT increased from north to south, likely due to longer periods of active evapotranspiration. A strong negative relationship appears to exist between Z and HT irresp ective of the soil class. Further this relationship appears to be stronger for SiL and SiC (Figure 2 5) compared to the remaining two classes. Field capacity also showed a negative relationship with HT for SL, SiL, and SCL though it was not as strong as in the case of Z , while for SiC there appeared to be no relation. A weak to moderate positive relationship between wp and HT appeared to be present, which increased as the percentage of clay increased from SL to SiC. All of the remaining parameters did not show any substantial relation to HT for any soil, except DC which showed a weak positive relationship with HT for S iC. Model output i.e. HT appears to vary linearly and monotonically for the soil textures under study. Sensitivity Indices Correlation a nalysis b ased s ensitivity i ndices Sensitivity indices based on partial correlation sensitivity analysis, namely PCCs and PRCCs, and the parameter rankings based on them are summarized in Table 2 3. The sign of PCC and PRCC indicate the nature of correlation between the two variables while the magnitude indicates the strength of the correlation. Root zone depth and fc showed negative sensitivity indices, indicating large values of these parameters are associated with small values of HT. On the other hand HT had positive sensitivity indices for DC, CN , and wp for all soils. Maximum uptake factor showed positive PCCs and PRCCs for SiL and SiC, while in the case of SL and SCL PCCs were negative while PRCCs were positive. However it should be noted that the corresponding magnitudes were relativ ely small. PAGE 55 55 The parameter importance rankings obtained from relative magnitudes, based on PCCs and PRCCs, are identical for SL, SiL, and SCL, where Z , fc , and wp are ranked first, second, and third, respectively, while DC , MUF , and CN came fourth, fifth, and sixth respectively. In the case of SiC parameter ranks differed from other soil classes with wp , and DC ranked second and third, respectively, followed by fc , CN , and MUF at four, five, and six respectively. ndices The first order and total effect sensitivity indices and parameter rankings (based on total effect indices) are presented in Table 2 4. Depending on soil type, Z contributed 63% to 92% to the total variability in HT and was ranked first for all soil types and location combinations. Note that for SiL and SiC the first order indices for Z (>90%) were higher than that for the other two classes, which was also seen from scatterplots. Soil hydraulic parameters, equivalent of fc and wp taken together, was the s econd most important parameter for three soil classes (SL, SiL, and SCL) at all locations and accounted for 4% to 35% of the variability in HT. In the case of SiC, SHP was ranked third at three locations where DC was ranked second. Amongst the remaining pa rameters, DC was ranked third, overall. Its contribution to the total variability was relatively small except for SiC where its contribution was similar to that of SHP . For Perrine, Saint Leo , and Gainesville (and SiC combination), DC was ranked second whi le SHP was ranked third. If the parameter rankings were to be based on first order effects, DC would have been ranked second at only one location for SiC. This is an indication of stronger interactive effects of DC compared to SHP for the case of SiC. The MUF and CN did not contribute substantially to the variability in HT ( S i and ST i < 1%) and they were ranked fourth and fifth, respectively, for SL, SiL, and SCL. In the case of SiC this order was interchanged with CN contributing more than 1%. This was not surprising as surface runoff generation would be expected to play a more important role in water balances for fine soils. PAGE 56 56 S i ), is an indicator of model type. For a perfectly additive S i = 1 or 100% method (Table 2 S i ranged from 91.4% to 99.6% indicating ARID is additive to a large degree, though not perfectly. Observations about model type can als o be made from PCCs, PRCCs, and the coefficient of model determination values on raw and ranked data (Table 2 5), by fitting a multilinear regression model. Saltelli et al. (2004) recommended R 2 raw 0.7 as a threshold to classify a model as linear and us ing PCCs as the sensitivity indices. Manache and Melching (2008) suggested looking at R 2 rank values as well. If they are greater than the corresponding R 2 raw , PRCCs are better indicators of sensitivity and the model is likely to be non linear and monotoni c. Table 2 5 shows that R 2 rank values are invariably greater than those of R 2 raw values, though all are greater than 0.7. The difference between R 2 raw and R 2 rank values vary from soil to soil with SL and SCL showing differences greater than those for SiL and SiC. Based on this it may be tenuous to say whether the model is linear or non linear. However, considering scatterplots (Figure 2 5) and the sum of first ord er sensitivity indices (Table 2 4) the model can be considered monotonic and linear. Effect of L ocation The locations used in the analysis had direct impact on HT as can be seen from the scatterplots, histograms , and other statistics. As noted earlier, HT increased from north to south. Location did not seem to impact PCCs or PRCCs substantially. Sensitivity indices obtained in rankings for Smithfield and Tifton we re the same. These ranks were different from those of Gainesville, Saint Leo , and Perrine, which showed the same parameter rankings. PAGE 57 57 Discussion ARID was originally developed for medium textured soils (sandy loam) and it is speculated that it may underesti mate the stress for sandy soils while overestimating for heavy soils when the default values of Woli et al. (201 2 ) are used. Si lty clay contains more than 40% clay and less than 20% sand and is the only truly heavy soil considered in this study. Though app ropriate parameter ranges were selected for this soil in this study , previous studies on ARID did not include heavy soils (e.g. SiC) in model evaluation, calibration, and validation studies. However inclusion or exclusion of SiC results from this analysis does not change the fact that , in general, Z was the most important parameter that influences ARID followed by soil hydraulic properties ( fc and wp or SHP ) irrespective of the method of analysis employed. In the case of SiC, DC has substantial effect on model output as well. Amongst the six paramters of ARID, two parameters ( Z and MUF ) were assigned the same distributions irrespective of soil type, while for DC the distibutions differed from soil to soil though the ranges were quite similar (Table 2 1) . The parameters fc , wp (or SHP ) , and CN were the only parameters which were assigned different probability distributions for different soils. Hence these parameters are likely to be responsible for the differences in the output distributions of different soil classes. However, CN which was assigned considerably different ranges for each soil type, seems to have very little influence over HT distributions (soil to soil variation). Thus we are left with fc and wp to explain this inter soil variability in HT . Statistics (minimum, maximum and central t endancy) of g enerated values of fc and wp (Table 2 6) did not relate well with HT distibutions . For example minimum valus of wp were similar for SL and SiL while maximum fc values of SiL and SCL were similar. However, CI limits and mean values of HT we re not similar for SL and SiL or SiL and SCL (Table 2 2). Also SiC, which had the highest values of wp and fc , had HT statistics somewhat similar to SCL. T o further explore PAGE 58 58 this soil type dependent variability in HT, values of total available water conte nt ( T AW C = fc wp ) were calculated (Table 2 6) and plotted against HT (Figure 2 6) . In comparing Figure 2 6 with the plots of wp and fc against HT (Figure 2 5), it can be seen that HT showed a stronger dependency on T AW C than to either wp and fc alone. This is intuitive because large T AW C means larger soil water storage capacity. In the work of Woli et al. (2012), the default value of T AW C was 0.13 which corresponds to medium textured soils (hydrologic group B and C soils) and was based on the st udies by Ratliff et al. (1983) and Ritchie et al. (1999). However, the soil databases used in those studies only partially cover these soil classes. Table 2 6 shows that for the two soils belonging to group B and C, T AW C mean values were considerably diffe rent than the default value previously used in ARID. Comparison W ith O ther S tudies Uncertainty and sensitivity analyses of ARID, other drought index models, and similar water balance approaches have been performed; e.g. Prats and PicÃ³ (2010) , Woli et al. (2013) etc. Typically these analyses were based on a one parameter at a time (OAT) approach, a type of LSA. Most of these studies considered a single location (different for different studies) and scenarios based on soil classes were generally not con sidered. In the studies in which different soil classes were considered, the impact of intra soil class variability was not accounted for. Five such studies are summarized in Table 2 7. A brief comparison of these studies with that of the current study is presented below. Woli et al. (2013) performed UA/SA of ARID using FAST by considering five locations (different from this study), four seasons, and three different PDFs for each of DC , CN , and TAWC . Average ARID was considered as the model output. A majo r difference between our study and Woli et al. (2013) is the consideration of soil classes. Woli et al. (2013) concluded that TAWC was the most important parameter followed by Z , MUF , DC , and CN in descending PAGE 59 59 order. TAWC and Z accounted for more than 80% of the total variability and the interactions were very small for all location season combinations. The order of ranking of these two parameters i.e. Z and TAWC (Woli et al., 2013) or SHP (this study) differed, as well as their contribution to variability differed substantially. Secondly, MUF which was ranked third in Woli et al. (2013) (contributing up to 15% to the total variability) was found to contribute less than 1% in our analyses and was ranked fourth or fifth. These differences in the results may b e the artifacts of considering soil classes in our analysis, consideration of different output, differences in estimated parameter distributions especially for Z , and to a minor extent the number of model runs used (257 used by Woli et al. (2013)). These t wo studies shall be looked at as complimentary rather than contradicting due to different experimental settings and the fact that the results match qualitatively with Z and TAWC or SHP being the important parameters. The uncertainty analysis of an irrigation scheduling model by Prats and PicÃ³ (2010) considered 12 textural soil classes and concluded that total available water ( TAW ), Z , and the fraction of T AW that can be depleted from the root zone before moisture stress be gins to occur ( p ) were the most important input parameters. Though they considered 12 soil textural classes, soil hydrualic properties were represented only by mean values for each soil class. Some other differences with the study by Prats and PicÃ³ (2010) were (a) they used different output criteria , though relative evaporation deficit , a term somewhat equivalent of ARID, was an integral part of the outputs, (b) they conducted simulations for a different crop ( o range) , and (c) UA/SA were perfo r med using an OAT approach considering only three parameters p , Z , and yield response factor to water stress (this parameter is not comparable to ARID parameters). Two UA/SA studies on AquaCrop, a canopy level crop model (Raes et al., 2009; Steduto et al., 2009) have previously been conducted. This model is water driven as the soil water PAGE 60 60 balance and water stress are at the c ore of all other processes. Being a crop mo del , it has many other crop specific parameters when compared to ARID. Salemi et al. (2011) applied this model in deficit irrigation management of winter wheat for a site in Iran. In their sensititvity analysis using a local, OAT method comparing simulated results to observed moisture contents, the model was found to be moderately sensitive to TAW and Z , while relatively less sensitive to initial water content. In a similar study of AquaCrop for q uinoa , Greets et al. (2009) found that the model was sensitiv e to Z , initial soil water content, fc and wp (apart from some other crop coefficients not related to water stress). As discussed earlier TAW is a product of TAWC and Z. Hence the importance of TAW and Z found by Prats and PicÃ³ (2010) and Salemi et al. (2011) potentially implies strong interactive effects between TAW and Z. However both studies employed an OAT method, and hence TAW should be considered as a proxy for TAWC . Padilla et al. (2011) developed a water balance model to simulate evapotranspiration for wheat and corn. However the model was more complex than ARID as it consider ed a number of crop coefficients. Further they performed sensitivity analysis by the OAT met hod considering only three parameters (two parameters related to root zone and one parameter similar to TAW ). A mongst these three parameters maximum rooting depth , which is equivalent to Z , was found to be the most important parameter . This results is sim ilar to that of the current study . In our study initial moisture content was not considered as an uncertain parameter though it was found to be an important parameter by Satti et al. (2004) and Greets et al. (2009) . Those studies were based on much shorte r time periods. Due to the length of the study period used in this analysis (24 years), initial soil moisture content was not likely to be an important parameter as it would only influence the very beginning of each model simulation. However, it may play a n PAGE 61 61 important role when smaller simulation periods are involved, for which its importance should be evaluated. Summary Global sensitivity and uncertainty analysis of ARID was successfully perfo r med in this study considering four textural soil classes and five locations across the southeast USA . One of the challenges was the presence of correlated model parameters when employing a variance based sensitivity technique. To account for this issue an innovative strategy, which comprised of sampling percent sand and clay from the soil parameter space to combine correlated parameters into a single parameter, was applied. Results of partial correlation analysis and the method of were found to be comparable . SA results indicate ARID is a monotonic and linear model. ARID was found to be most sensitive to Z (ranked first in both analyses) followed by soil hydraulic properties i.e. SHP or fc and wp indicating that these properties are important irrespective of the soil class. These results are consistent with with the conclusions of previous studies on ARID and soil water balance models. Extra care may be required while using this model for clayey so ils since DC was found to increase in importance for heavier soils. In many crop and drought index models and applications, soil hydraulic properties are deterministically parameterized based on their soil class. The variability of ARID i.e. uncertainty in outputs mainly depends on soil type (smaller ranges for finer soils against wider ranges for soils with greater sand content) . Hence it will be advisible to determine soil hydraulic properties in detail instead of specifying soil type alone whereever p ossible, particularly in case of coarser soils. Weather location did not have substantial influence on parameter importance rankings except for SiC. PAGE 62 62 This analysis covered all four soil hydrologic groups. However, only one soil class belonging to each gro up was considered. Hence, directly extending the results of this analysis to other soil classes without verification should be done with caution. PAGE 63 63 Table 2 1 . List of input parameters and summary of their estimated/ proposed probability distribution funct ions Parameter 1, 2 Units Range Distribution Distribution Parameters 3 Z mm 250 800 Uniform a=250, b=800 MUF day 1 0.075 0.11 Normal Sandy Loam DC day 1 0.35 0.65 Triangular a=0.2, b=0.6, c=0.7 CN 30 50 Uniform a=30, b=50 wp m 3 m 3 (%) 5.5 13 Uniform a=5.5, b=13 fc m 3 m 3 (%) 11 27 Normal SHP Discrete Uniform n=100 Silty Loam DC day 1 0.2 0.6 Uniform a=0.2, b=0.6 CN 50 67 Uniform a=50, b=67 wp m 3 m 3 (%) 5.5 17.3 Uniform a=5.5,b=17.3 fc m 3 m 3 (%) 20.5 36.9 Triangular a=20.5, b=34, c=36.9 SHP Discrete Uniform n=100 Sandy Clay Loam DC day 1 0.25 0.65 Uniform a=0.2, b=0.7 CN 67 77 Uniform a=67, b=77 wp m 3 m 3 (%) 14 22 Uniform a=14, b=22 fc m 3 m 3 (%) 21.2 33.8 Beta SHP Discrete Uniform n=100 Silty Clay DC day 1 0.2 0.4 Uniform a=0.1, b=0.45 CN 77 85 Uniform a=77, b=85 wp m 3 m 3 (%) 24 31.9 Triangular a=24, b=24.5, c=31.9 fc m 3 m 3 (%) 38.6 43.6 Triangular a=38.6, b=40, c=43.6 SHP Discrete Uniform n=100 1 Z root zone depth, MUF water uptake factor, DC drainage coefficient, CN curve number, wp wilting point, fc field capacity, SHP soil hydraulic properties 2 wp and fc distributions are used for correlation based SA while SHP is used for Sobo 3 for triangular distribution: a = minimum, b = typical value and c = maximum value , Uniform distribution: a = minimum, b = maximum , , Beta , Di screte Uniform: n = number of discrete points PAGE 64 64 Table 2 2 . Characteristics of the distributions of high threat (HT) 1 (LL=lower limit, UL=upper limit) for the four soils used in the study (SL= Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) Location SL SiL SCL SiC Mean LL UL Mean LL UL Mean LL UL Mean LL UL Smithfiled, NC 0.12 0.06 0.27 0.06 0.02 0.11 0.11 0.06 0.22 0.08 0.05 0.14 Gainesville, GA 0.15 0.08 0.28 0.08 0.04 0.14 0.14 0.09 0.23 0.11 0.07 0.16 Tifton, GA 0.20 0.12 0.36 0.12 0.07 0.18 0.19 0.13 0.30 0.15 0.11 0.22 Saint Leo, FL 0.22 0.13 0.39 0.13 0.08 0.21 0.21 0.14 0.33 0.17 0.12 0.25 Perrine, FL 0.23 0.14 0.38 0.14 0.08 0.22 0.22 0.15 0.33 0.18 0.13 0.25 1 HT corresponds to ARID values > 0.6 PAGE 65 65 Table 2 3 . Partial Correlation Coefficient (PCC), Partial R ank Correlation Coefficient (PR C) and Parameter Rankings Smithfield, NC Gainesville, GA Tifton, GA Saint Leo, FL Perrine, FL Sandy Loam Para 2 SI 1 Ranks SI Ranks SI Ranks SI Ranks SI Ranks PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C DC 0.24 0.35 4 4 0.28 0.41 4 4 0.28 0.37 4 4 0.32 0.39 4 4 0.36 0.41 4 4 Z 0.90 0.97 1 1 0.91 0.97 1 1 0.92 0.97 1 1 0.93 0.97 1 1 0.94 0.97 1 1 CN 0.01 0.00 6 6 0.01 0.00 6 6 0.01 0.00 6 6 0.01 0.00 6 6 0.01 0.00 6 6 MU F 0.05 0.10 5 5 0.04 0.13 5 5 0.03 0.12 5 5 0.06 0.04 5 5 0.06 0.03 5 5 wp 0.69 0.87 3 3 0.71 0.86 3 3 0.73 0.87 3 3 0.77 0.87 3 3 0.79 0.87 3 3 fc 0.75 0.90 2 2 0.77 0.90 2 2 0.79 0.90 2 2 0.82 0.90 2 2 0.84 0.90 2 2 Silty Loam SI Ranks SI Ranks SI Ranks SI Ranks SI Ranks PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C DC 0.31 0.35 4 4 0.45 0.43 4 4 0.40 0.38 4 4 0.42 0.42 4 4 0.52 0.47 4 4 Z 0.97 0.99 1 1 0.98 0.99 1 1 0.98 0.99 1 1 0.98 0.99 1 1 0.99 0.99 1 1 CN 0.03 0.04 6 6 0.05 0.05 6 6 0.06 0.06 6 6 0.07 0.08 6 6 0.06 0.06 6 6 MU F 0.16 0.24 5 5 0.26 0.32 5 5 0.37 0.37 5 5 0.11 0.16 5 5 0.07 0.09 5 5 wp 0.83 0.89 3 3 0.88 0.89 3 3 0.88 0.89 3 3 0.87 0.89 3 3 0.91 0.89 3 3 fc 0.84 0.90 2 2 0.89 0.90 2 2 0.89 0.90 2 2 0.88 0.90 2 2 0.92 0.90 2 2 Sandy Clay Loam SI Ranks SI Ranks SI Ranks SI Ranks SI Ranks PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C DC 0.32 0.52 4 4 0.37 0.58 4 4 0.37 0.55 4 4 0.42 0.57 4 4 0.46 0.60 4 4 Z 0.88 0.98 1 1 0.89 0.98 1 1 0.90 0.98 1 1 0.92 0.98 1 1 0.93 0.98 1 1 CN 0.04 0.08 6 6 0.04 0.08 6 6 0.05 0.11 5 6 0.06 0.12 5 5 0.05 0.09 6 5 MU F 0.05 0.13 5 5 0.04 0.16 5 5 0.03 0.17 6 5 0.06 0.06 6 6 0.06 0.05 5 6 wp 0.70 0.90 3 3 0.71 0.89 3 3 0.73 0.89 3 3 0.76 0.90 3 3 0.78 0.90 3 3 fc 0.77 0.93 2 2 0.78 0.93 2 2 0.80 0.93 2 2 0.82 0.93 2 2 0.84 0.93 2 2 Silty Clay SI Ranks SI Ranks SI Ranks SI Ranks SI Ranks PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C PCC PR C DC 0.55 0.65 3 4 0.60 0.71 3 3 0.60 0.66 3 3 0.65 0.73 3 3 0.74 0.78 3 3 Z 0.97 0.99 1 1 0.97 0.99 1 1 0.97 0.99 1 1 0.97 0.99 1 1 0.98 0.99 1 1 CN 0.45 0.63 5 5 0.44 0.62 5 5 0.50 0.63 5 5 0.50 0.64 5 5 0.45 0.52 5 5 MU F 0.22 0.40 6 6 0.25 0.44 6 6 0.30 0.46 6 6 0.10 0.18 6 6 0.09 0.14 6 6 wp 0.74 0.84 2 2 0.72 0.83 2 2 0.76 0.82 2 2 0.76 0.84 2 2 0.82 0.85 2 2 fc 0.53 0.67 4 3 0.51 0.65 4 4 0.55 0.64 4 4 0.56 0.67 4 4 0.63 0.68 4 4 1 SI Sensitivity Index 2 DC drainage coefficient (day 1 ), Z root zone depth (mm), CN curve number, MUF max. uptake factor (day 1 ), wp soil water content at wilting point (m 3 m 3 ), fc soil water content at field capacity (m 3 m 3 ) PAGE 66 66 Table 2 4 . Smithfield, NC Gainesville, GA Tifton, GA Saint Leo, FL Perrine, FL Sandy Loam Parameter 1 S i 2 ST i 3 Rank S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank DC 0.9% 2.7% 3 1.2% 3.0% 3 1.1% 2.7% 3 1.2% 2.8% 3 1.3% 2.9% 3 Z 63.3% 71.3% 1 62.9% 70.8% 1 64.0% 70.3% 1 65.3% 70.1% 1 65.3% 69.6% 1 CN 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 MUF 0.1% 0.5% 4 0.1% 0.6% 4 0.1% 0.5% 4 0.1% 0.4% 4 0.1% 0.3% 4 SHP 27.1% 34.9% 2 27.7% 35.3% 2 28.0% 34.2% 2 28.4% 33.1% 2 28.9% 33.2% 2 Silty Loam S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank DC 0.4% 0.4% 3 0.7% 0.7% 3 0.3% 0.4% 3 0.5% 0.6% 3 0.8% 0.9% 3 Z 90.3% 91.6% 1 89.2% 90.0% 1 89.4% 90.4% 1 89.2% 90.3% 1 88.0% 89.8% 1 CN 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 0.0% 0.0% 5 MUF 0.0% 0.2% 4 0.2% 0.3% 4 0.4% 0.4% 4 0.0% 0.1% 4 0.0% 0.0% 4 SHP 8.9% 8.9% 2 8.9% 9.5% 2 9.0% 9.1% 2 9.1% 9.2% 2 9.7% 9.8% 2 Sandy Clay Loam S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank DC 2.2% 3.3% 3 2.7% 4.1% 3 2.4% 3.5% 3 2.6% 3.7% 3 3.1% 3.9% 3 Z 82.8% 85.8% 1 82.1% 84.8% 1 82.8% 85.1% 1 83.7% 85.3% 1 83.5% 84.8% 1 CN 0.0% 0.1% 5 0.0% 0.0% 5 0.0% 0.1% 5 0.0% 0.1% 5 0.0% 0.0% 5 MUF 0.1% 0.5% 4 0.1% 0.7% 4 0.1% 0.5% 4 0.2% 0.4% 4 0.1% 0.3% 4 SHP 10.8% 13.0% 2 10.8% 12.9% 2 10.9% 12.5% 2 11.0% 12.2% 2 11.3% 12.2% 2 Silty Clay S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank S i ST i Rank DC 2.4% 3.7% 3 3.2% 4.7% 2 2.7% 3.8% 3 3.4% 4.6% 2 4.1% 5.0% 2 Z 91.3% 92.5% 1 89.8% 91.2% 1 90.1% 91.3% 1 90.5% 91.5% 1 90.8% 91.5% 1 CN 1.4% 1.5% 4 1.4% 1.5% 4 1.7% 1.8% 4 1.6% 1.6% 4 0.8% 0.9% 4 MUF 0.2% 0.3% 5 0.4% 0.5% 5 0.5% 0.6% 5 0.0% 0.0% 5 0.0% 0.1% 5 SHP 3.8% 4.5% 2 3.6% 4.1% 3 3.9% 4.3% 2 3.7% 4.1% 3 3.9% 4.2% 3 1 DC drainage coefficient (day 1 ), Z root zone depth (mm), CN curve number, MUF max. uptake factor (day 1 ), SHP soil hydraulic properties 2 S i First order sensitivity index 3 ST i Total effect sensitivity index PAGE 67 67 Table 2 5 . Coefficient of model determination values (R 2 ) based on raw and ranked output obtained in partial correlation analysis ( for the four soils used in the study (SL = Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) SL SiL SCL SiC R 2 Raw R 2 Rank R 2 Raw R 2 Rank R 2 Raw R 2 Rank R 2 Raw R 2 Rank Smithfield, NC 0.8498 0.9589 0.9512 0.9753 0.8402 0.9656 0.9459 0.9769 Gainesville, GA 0.8618 0.9577 0.9667 0.9749 0.8492 0.9644 0.9396 0.9748 Tifton, GA 0.8770 0.9584 0.9681 0.9745 0.8651 0.9648 0.9504 0.9733 Saint Leo, FL 0.9015 0.9590 0.9649 0.9753 0.8872 0.9656 0.9530 0.9769 Perrine, FL 0.9139 0.9589 0.9751 0.9752 0.9008 0.9654 0.9671 0.9785 PAGE 68 68 Table 2 6 . Summary of wilting point ( wp ), field capacity ( fc ), and Total available water content ( TAWC ) (SL= Sandy Loam, SiL = Silty Loam, SCL = Sandy Clay Loam and SiC = Silty Clay) wp fc TAWC Min Max Mean 1 Min Max Mean Min Max Mean SL 0.06 0.13 0.09 0.11 0.27 0.19 0.05 0.16 0.10 SiL 0.06 0.17 0.11 0.21 0.37 0.30 0.11 0.27 0.19 SCL 0.14 0.22 0.18 0.21 0.34 0.28 0.04 0.16 0.10 SiC 0.24 0.31 0.27 0.39 0.44 0.41 0.09 0.17 0.14 1 Median values of wp fc and TAWC are similar to mean values PAGE 69 69 Table 2 7 . Summary of uncertainty and sensitivity analyses of ARID and other soil water balance models Source Model Method Location Meteorological Data 1 Soil Classes 2 Important Parameters 3 , 4 Woli et al. (2013) ARID GSA FAST Multiple No No (1) TAWC (2) Z (3) MUF and (4) DC Prats and PicÃ³ (2010) Irrigation scheduling soil water balance approach based on FAO 56 LSA OAT Single Yes Yes 5 TAW, Z, p Greets et al. (2009) AquaCrop LSA OAT Single Yes No fc wp , initial soil water content Salemi et al. (2011) AquaCrop LSA OAT Single Yes No TAW, Z Padilla et al. (2011) A soil water balance model LSA OAT No No Z max 1 Whether meteorological data was treated as probabilistic: Yes or No 2 Whether soil classes were considered: Yes or No 3 TAWC total available water content (m 3 m 3 ), Z root zone depth (m), MUF maximum uptake factor (day 1 ), DC drainage coefficient (day 1 ), TAW total available (m), p fraction of TAW that can be depleted from the root zone before moisture stress begins, fc soil water content at field capacity ( m 3 m 3 ), wp soil water content at wilting point (m 3 m 3 ), Z max maximum rooting depth (m) 4 Only parameters similar to those of ARID were considered 5 Though soil classes were considered intra soil class variability was not considered PAGE 70 70 Figure 2 1 . Schematic Soil Water Balance . The square represents the maximum volume of water stored in the soil, with three moisture regions from easily extractable to more difficult (saturation, field capacity, wilting point) (adopted from Allen et al. 1998) PAGE 71 71 Figure 2 2 . Locations of weather stations used in this study PAGE 72 72 Figure 2 3 . Long term mean monthly precipitation, maximum temperature, and minimum temperature at A ) Smithfield, NC, B ) Gainesville, GA, C ) Tifton, GA, D ) Saint Leo, FL, and E ) Perrine, FL. (based on 1980 2003 daily data) A B C D E PAGE 73 73 Figure 2 4 . Histograms of HT method PAGE 74 74 A SL B SiL PAGE 75 75 Figure 2 5 . Scatterplots of HT against model parameters A ) Sandy Loam (SL), B ) Silty Loam (SiL), C ) Sandy Clay Loam (SCL), and D ) Silty Clay (SiC). C SCL D SiC PAGE 76 76 Figure 2 6 . Scatterplots of HT against Total Available Water Content ( TAWC ) for all soil types and location combinations PAGE 77 77 CHAPTER 3 A MULTI CRITERIA TRAJEC TORY BASED PARAMETER SAMP LING STRATEGY FOR THE SCREENING METHOD OF ELEMENTARY EFFECT S Environmental models are increasingly used in a regulatory decision making framework to tackle a variety of complex problems (e.g. Jakeman et al., 2006; Warmink et al., 2 010 ; Bennett et al., 2012 ). Both the European Commission (EC, 2009) and the United States Environmental Protection Agency (EPA, 2009) have emphasized the importance of risk assessment and have provided the guidelines on assessing and reducing the risks in environmental model applications. Sensitivity analysis is one of the fundamental tools used for model evaluation (Saltelli et al., 2000; Refsgaard et al., 2007; Liu et al., 2008; Saltelli et al., 2008). Today, several sensitivity analysis methods are available ranging from derivative based local approaches to more rigorous variance based global sensitivity analysis (GSA) methods such as the Fourier Amplitude Sensitivity Test (FAST) and the Cukier et al., 1978 ; Sobol , 1993 ; Saltelli e t al., 2000; Helton and Davis, 2002; Cacuci and Ioneschu Bujor, 2004; Saltelli et al. , 200 5 ). The merits of GSA methods over local methods are well established (Saltelli et al., 2004; Saltelli et al., 2008) . However, variance based GSA methods become time consuming with an increase in the number of parameters to the point of computational infeasibility ( Saltelli et al., 2005 ; Kucherenko et al., 2009 ) . This is often the case with many environmental and ecological models as they are typically characterized by tens to hundreds of parameters where the consideration of spatial variability further increases the model dimensionality, ultimately requiring a large number of model runs and computational resources ( van Griensven et al., 2006; Jawitz et al., 2008; Fogli a et al., 2009; Makler Pick et al., 2011; Ciric et al., 2012; Herman et al., 2013a,b; Moreau et al., 2013 ). Thus parameter screening, the initial separation of important model parameters from unimportant ones, becomes essential before applying rigorous var iance based GSA methods. PAGE 78 78 Various methods such as first and higher order derivatives, one at a time (OAT) approaches , the method of elementary effects (EEs) (Morris, 1991), the systematic fractional replicate design (Cotter, 1979), the iterated fractional f actorial design (Andres and Hajas, 1993) and the sequential bifurcation design (Bettonvil, 1990; Bettonvil and Kleijnen, 1997) have been developed for low computational cost parameter screening exercises . Among these, the method of EEs (and its variants) i s the most widely used screening method in environmental modelling studies and has been recommended as a part of a modern statistical framework for GSA (e.g. Campolongo et al., 1999; Cariboni et al., 2007; Carpena et al., 2007 ; Tong and Graziani, 2007; Yang et al., 2012; Zhan et al., 2013 ). However, the original method of EEs i.e. Morris (1991) method was found to have drawbacks regarding parameter sampling and the calculation of sensitivity measures sometimes leading to unreliable parameter screen ing (Campolongo et al., 2007). Santer et al. (2003) categorized parameter sampling schemes for computer experiments into three types: (1) Latin Hypercube Sampling ( LHS ) , (2) distance measure criteria , and (3) uniformity . Each type has advantages and disad vantages (e.g. Fang et al., 2000; Johnson et al., 1990; Santiago et al., 2012a,b; Wiens, 1991 ). However, a sampling strategy belonging to only one of the above three types is not totally satisfactory for sampling designs (Santer et al., 2003). They also su ggested that care should be taken when combining design types to avoid computational infeasibility. It should be noted that most of the computer experiment sampling strategies generate parameter samples in a unit hyperspace which are then transformed to us er defined distributions. Table 3 1 summarizes the original and recent developments in parameter sampling and sensitivity measures for the EEs method. Theoretical basis of these strategies vary considerably. PAGE 79 79 Each successive work has sought to achieve a mor e effective sampling strategy, i.e. exploration implies the selection of fewer trajectories for reliable screening as well as the time required for sample generatio n. Morris (1991) discussed advantages of trajectory based designs over LHS when one needs to restrict the number of model runs. The stratified sampling used by Morris (1991) from a definite number of levels ( p ) did not employ a distance measure or distribu tion criteria and can be regarded as a random sampling over the gridded parameter space . This sampling strategy by Morris will be referred as MM hereon in this article. While LHS may not be eff icient, some of its variants have attractive properties such as orthogonality , which can be useful for screening methods (Santer et al., 2003). Van Griensven et al. (2006) presented a LHS based sampling for the EE method while a simplex based sampling strategy was introduced by Pujol (2009). The strategy of optimized trajectories (OT) (Campolongo et al., 2007) and the modified optimized trajectories (MOT) (Ruano et al., 2012) base d their sampling on maximizing Euclidean distance ( ED ) between trajectories. Note that although ED is not the only possible distance criteria, it is the only one explored so far for sample optimization ( Santer, 2003 ; Campolongo et al., 2007 ). The c ell based trajectory sampling (Saltelli et al., 2009) uses a rather complicated parameter sampling scheme. Given an individual input parameter X i , EEs are calculated by tracking and combining model values for sample points which differ from each other in all parameter coordinates except the X i th . This method, originally designed for models with strong interact ive effects, was found to be inferior to some other methods in identifying important parameters based on total sensitivity (Campolongo et al., 2011). The radial sampling strategy proposed by Campolongo et al. (2011) was designed as an extension of the basi c OAT PAGE 80 80 quasi rameters is high (Saltelli et al., 2000) and repeated use of the trajectory base points makes success of this strategy dependent on specific points (see Herman et al., 2013a). It is noteworthy that n one of the screening sampling strategies discussed above were based on the principle of reproducing uniform parameter distribution s closely . A uniform parameter distribution is indicated in this context by the probability density function of an individual parameter values from all trajectories. However, only tw o of the above strategie s (Campolongo et al., 2007 and Ruano et al., 2012) discuss the importance of matching the distributions of generated parameter samples. One of the possible reasons that uniform distribution sampling has not been previously considere d could be the associated run time efficiency, as mentioned in Santer el al. (2003). In addition, n one of the previous studies , with the notable exception of Ruano et al. ( 2012) , consider ed the issue of run time , which is an important feature from the application point of view. The MOT method was developed by Ruano et al. (2012) in response to the large computational time requirement for sample generation using the OT, especially for high dimensional models. In spite of the availability of an array of sampling strategies for EE based parameter screening , a n efficient sampling strategy which is suitable across multiple criteria for all types of models still eludes the modeling community. The objective of this research was to develop a new screening sampl ing strategy by combining t wo parameter sampling criteria , so that it (i) can be used for a wide range of models, (ii) is fast enough from a practical point of view , and (iii) produces parameter distributions that closely resemble the theoretical ones . The performance of this newly developed sampling strategy, the strategy of sampling for uniformity (SU), was PAGE 81 81 compared to those of MM , OT, and MOT popular trajectory based sampling strategies. The various evaluation criteria include d (i) the quality of the g enerated parameter distribution s , (ii) ED of the sampling trajectories , (i ii ) time for sample generation, and ( i v) the ability to correctly identify important model parameters i.e. parameter screening efficiency . Material and Methods Method of EEs The me thod of EEs is based on the concept of building trajectories in a unit hypercube parameter space . C oordinates of two successive points within a trajectory differ from each other only in one dimension , or one parameter coordinate , by a fix ed amount . Note that f or the case of radial sampling strategy is variable and coordinates of all sample points in a trajectory differ from the base point in exactly one dimension . is a multiple of 1/( p 1), usually set to p /[2( p 1)], where p is the number of levels alo ng any parameter axis. p levels indicate that each parameter can take values from the set {0, 1/( p 1), 2/( p ure 3 1 A ). A common practice is to use four levels ( p = 4 ; e.g. Campolongo et al., 2007; Saltelli et al., 2008; Ruano et al., 2012 ). S ome sensitivity analysis tools, e.g. SimLab v2.2 ( Joint Research Center, http://ipsc.jrc.ec.europa.eu/?id=756 ), use truncated unit hyperspace (Fig ure 3 1 B ) in order to avoid infeasibilities in converting trajectories from unit hyperspace to real values for probability distributions with infinite tails, where is scaled accordingly. The EE associated with the i th parameter is calculated using (3 1) where y represents the model under study and k is the number of model parameters. One EE is produced per parameter from each trajectory. Since the sample consists of r trajectories, r EEs for each parameter are obtained from the entire sample. Morris (1991) used the mean and PAGE 82 82 standard deviation of the EEs of the i th parameter as the global sensitivity measures. However, for non monotonic models the method of EEs may not detect some parameters to be influential due to positive and negative EE values canceling ea ch other to produce negligible . To address this, Campolongo et al. (2007) recommended calculating the mean of the absolute values of the EEs ( * ). The values of * (or ) and have been found equivalent to total effect and interaction effect sensitivity indices obtained in variance decomposition sensitivity analyses ( Morris, 1991; Campolongo and Saltelli , 1997; Campolongo et al., 2007 ). Trajectories Based on Sampling for Uniformity (SU) The SU strategy proposed in this work is based on the principle of generating discrete uniform distributions (based on the first and the last points of trajectories) for all parameters and maximizing the spread of the sample. It is assumed that p = 4 and r is even, a common practice in EE screening exercises (Saltelli et al., 2008). The algorithm to generate a parameter sample for k parameters consisting of r trajectories, each trajectory containing k+ 1 points, by maximizing the ED between them from Q repe titions, is explained below with an example for the case r = 4, k = 2, p = 4, and Q = 300 (Fig ure 3 2 A to C ). Appendix A (Figure A 1) presents the Matlab code developed for SU which is available as a free download from ( http://abe.ufl.edu/carpena/software/SUMorris.shtml ). STEP 1: Since we desire to generate discrete uniform distributions for all parameters, the 1 st point of exactly r /4 or ( r /2Â±1)/2 trajectories should have a value on each parameter level, when r /2 is even or r /2 is odd, respectively. First we randomly select trajectories, for each parameter, for which the 1 st point falls on the first level. From the remaining trajectories, r / 4 or ( r /2Â±1)/2 trajectories are selected randomly for which the 1 st points have a value on the second level. The same procedure is repeated for the third level. The remaining trajectories at this stage take on the fourth level value at their 1 st point. Con sidering the property that the 1 st and the last PAGE 83 83 i.e. ( k +1) th points of a given trajectory are formed by levels with constant offset (Fig 3 1 A , B ), the last points of all trajectories are assigned appropriate parameter values. In the case of any duplicati on, the above procedure is repeated until unique points are sampled. At the end of this step coordinates of the 1 st and the ( k +1) th points of all trajectories are generated (Fig ure 3 2 A ). STEP 2: Intermediate points i.e. 2 nd to k th points of each trajecto ry are formed by changing one parameter coordinate of the previous point in that trajectory. This is done by randomly generating a perturbation vector of length ( k 1) for each trajectory, whose elements are non repeated integers from 1 to k . If the first e lement of this vector is 1 then only the first coordinate of the 1 st point of that trajectory is changed by to form the 2 nd point (Fig ure 3 2 B ). A complete parameter sample is generated at the end of this step. A check is provided for the uniqueness of intermediate points and in the case of duplication the above procedure is repeated till all sample points are unique (Fig ure 3 2 B ). STEP 3: Steps 1 and 2 are repeated Q times. The ED , as defined in equation 3 2 and 3 3 , between trajectories is calculated a t the end of each repetition. The trajectory set corresponding to the repetition for which ED = E D max is selected as the sample for EE analysis (Fig ure 3 2 C ). It should be noted that the distance maximization approach used in this procedure is more simple than that used for OT and MOT. Numerical Experiments Numerical experiments were carried out to compare performances of the four sampling strategies; namely MM, OT, MOT and SU, based on the four criteria; fitness to discrete uniform distributions, trajector y spread ED , computational efficiency CPU time requirement, and parameter screening efficiency. For the last criterion five standard test functions (Table 3 2 and Appendix B) with varying degrees of complexity and a range of k were used. T est functions include d: (a) the , (b) the modified G or G* function (Saltelli et PAGE 84 84 al., 2010) , (c) the Morris function or M function (Morris, 1991) , (d) the O function (Oakley and , and (e) the B function (Saltelli et al., 20 08) . These functions were selected based on (1) application history, (2) their nature (properties like linearity, monotonicity, curvature etc.), (3) the availability of analytical expressions for (or known values of) the total order sensitivity indices, an d (4) their modifiability to any k (except for O and M). The B function was used in four configurations while the G and G* functions were used in six configurations (Table 3 2). The M and O functions, by definition, can have only 20 and 15 parameters, resp ectively. Note that for all test functions parameters were uniformly distributed within 0 1. Thus, the complete set of numerical experiments consisted of combinations of the following characteristics p = 4, r = {6, 10}, k = {15, 20, 35, 50, 80, 100}, N = 5 00, Q = {1, 300, 1000} Where N is the number of times each experiment is repeated. By definition Q = r for MM , while for OT and MOT Q = 1000 is a recommended oversampling value (Campolongo et al., 2007; Ruano et al., 2012). For SU the ideal value of Q is u nknown and hence three Q values ( Q = 1, 300 and 1000), hereafter denoted as SU1, SU300 and SU1000, were used. For this study truncated parameter distributions (Fig ure 3 1b) similar to the MM implementation in SimLab v2.2 (JRC, http://ipsc.jrc.ec.europa.eu/?id=756 ) were used. Note that the meaning of Q (i.e. oversampling size) is different for SU and other sampling strategies. While in the case of MM, OT, and MOT Q indicates the size of the trajectory pool from which the final trajectories are selected, for SU Q indicates the number of times the whole sample generation procedure is repeated. Testing for uniformity of the generated distributions All sampling strategies for the EE method assume the parameters have discrete uniform distributions (in unit hyperspace). To check how consistent the generated distributions were, the PAGE 85 85 distributions were tested considering only the first and the last points of trajectories (e.g. C ampolongo et al., 2007) against the discrete uniform distribution using the Chi square goodness of fit test at three significance levels ( = 0.2, 0.1 and 0.05) . The f raction of model parameters failin g the Chi square goodness of fit test was reported. Tim e comparison Simulation time i.e. CPU time required to generate the final sample consisting of r trajectories , in seconds, required by all six sampling strategies were recorded and plotted against k for both r = 6 and 10 . All numerical experiments, except those for time requirements, were carried out using the High Performance Computing (HPC) Center Facilities at the University of Florida ( http://www.hpc.ufl.edu/ ) . Time requirement experiments were carried out on a de sktop machine with an Intel Core 2 Duo (3 GHz) processor with 8 GB RAM in order to avoid run time variations likely to occur on the HPC grid. Euclidean distance a b ( 3 2) where X i a ( z ) and X j b ( z ) are the z th coordinates of the i th and j th a b respectively. The Euclidean distance for the entire parameter sample is calculated as ( 3 3) Distributions of E D for combinations of r and k were compared by visual inspection and non parametric statistical tests. The Kruskal Wallis and the Wilcoxon rank sum tests (Wilcoxon, 1945 ; Kruskal and Wallis, 1952; Helsel and Hirsch, 2002 ) were used at the 5% significance level to test the separation o f medians. PAGE 86 86 Parameter screening efficiency Effectiveness of sampling strategies was assessed by calculating two skill scores, g (Campolongo et al., 2011) and R (defined later in this article) which were then integrated across all functions using two approac hes. Both the skill scores are different forms of the ratio representing the ability of EE screening ( * ) to correctly identify important model parameters (Eqs. 3 4 and 3 5 ). The g score merely represents the degree of match between sets of important param eters identified by screening and benchmark methods while the R score takes into account the one to one correspondence of ranks of the top parameters. Hence, R score is much more stringent indicator of screening efficiency if comparison is to be made betwe en the EE and variance based GSA methods. Parameters with total effect sensitivity index ( ST i ) > 5%, based on analytical expressions or prior knowledge, were assumed to be true important. ( This criterion was not strictly followed e specially when the important parameters were small in number. Please refer to Appendix B for the list of parameters identified as important for each function configuration used in this study ). For each function configuration the set of true important param eters ( T ST ), number of parameters in T ST ( #T ), and their ranks (Table 3 2 and Appendix B) were recorded. From screening experiments top #T parameters based on * (i.e. set T ) were recorded. Average g scores (avg_ g ) and average R scores (avg_ R ) for functi on configurations, r , and sampling strategy combinations were determined from the N = 500 runs. ( 3 4 ) ( 3 5 ) PAGE 87 87 For example, if T ST and T , with parameters in the descending order of importance are { X 1, X 2, X 3, X 4} and { X 1, X 3, X 2, X 4}, respectively, then the skill scores are: g = 4/4 = 1, R = 2/4 = 0.5. Approach 1: Average skill scores for a given function type (Avg_ g and Avg_ R ) (averaged over all configurations of a function type and r ) for each sampling strategy were calculated first. These ensemble averages were then further averaged across all function types to obtain A VG_ g and AVG_ R skill scores for sampling strategies. Ple ase refer equations 3 6 to 3 11 . Approach 2: This approach is similar to approach 1 but is founded on ranking the sampling strategies based on avg_ g and avg_ R scores. For each function configuration r combination sampling strategy having the highest avg_ g and avg_ R were ranked 1 for respective scores and the ones with the lowest were ranked 6. These ranks were denoted as avg Rank_ g and avg Rank_ R . These were then averaged for each function type to calculate Avg Rank_ g and Avg Rank_ R . Averaged ranks were fu rther averaged across all function types for each sampling strategy to calculate AVG R ank_ g and AVG R ank_ R scores. T he sampling strategy with the smallest AVG Rank score is the best and vice versa. Please refer eqs 3 12 to 3 17 . Note that Approach 1 may produce biased results if a sampling strategy consistently performs poorly or exceptionally well for one or more functions, which happens to be the case in this study as discussed later. Approach 2 can overcome this issue but involves ranking sampling stra tegies as a relative measure whereas Approach 1 which is an absolute measure. Equations for g score ( 3 6 ) PAGE 88 88 ( 3 7 ) ( 3 8 ) Equations for R score ( 3 9 ) ( 3 10 ) ( 3 11 ) Equations for g score ( 3 12 ) ( 3 13 ) ( 3 14 ) Equations for R score ( 3 15 ) ( 3 16 ) ( 3 17 ) w here r = the number of trajectories, k = the number of parameters, FT = the function type, FT k = the function configuration with k parameters, SS = the sampling strategy, N = the number of times screening experiment is repeated for each r , FT k , SS combination ( N = 500 for this study), PAGE 89 89 n 1 = the number of function configuration r combinations for the function type FT ( n 1 = 12, 12, 8, 2, 2 for G, G*, B, M and O, respectively), n 2 = number of function types ( n 2 = 5 for this study) R esults Uniformity of Generated Parameter Distributions The average fraction of parameters failing the Chi Square test for e ach sampling strategy is presented in Table 3 3. For a given k , MM had the highest failure rate followed by MOT and OT, respectively, while SU (all three configurations) had no failures. Irrespective of the sampling strategy (except SU) failures increased with k indicating performance deterioration as model dimensionality increased. SU having no failures was expected as the method was specifically designed to generate discrete uniform distributions (considering the first and the last points of sample trajec tories). Sampling Time Requirement The time required by the sampling strategies is plotted as a function of k in Figure 3 3. For MM, OT, and MOT as k increased from 15 to 100 the time requirement increased by factors ~ 1500, ~ 80 and ~ 200, respectively. In the case of SU, as k increased from 15 to 100 the increase in sampling time was merely ~30 40 times indicating the marked improvement in sampling e fficiency of this strategy in high dimensional parameter space. For all sampling strategies the time requirements increased non linearly with k . For r = 6 and 10 OT required the maximum time. On the other hand, SU1 took the least time for all k > 20 at r = 6. At r = 10 SU1 was the most time efficient only when k > ~ 27, otherwise MM was the most efficient strategy. SU300 required approximately two orders of magnitude greater time when compared to MM for k < 80. For higher k the difference was merely an orde r of magnitude. On the other hand OT and MOT required an order of magnitude time longer than SU300 (except MOT at k < 35). SU1000 needed similar times as those of MOT and OT, slightly higher than MOT but lower than OT, for PAGE 90 90 k = < 20. For higher k (> 35) SU1 000 was more efficient than OT and MOT. For MM, OT and MOT the time requirements did not increase much as r increased from 6 to 10 while for SU (all three configurations) the time requirement increased by ~ 2, which is approximately proportional to the inc rease in r (10/6). In terms of absolute time requirements, SU1 and MM needed just a few seconds while SU300 took a few minutes irrespective of k . For OT and MOT the time requirement varied from minutes ( k < 50) to a few hours ( k > 50). SU1000 required time on the order of hours only for k > 80. Spread of Trajectories/ Euclidean Distances A summary of the spread of trajectories i.e. ED in the parameter hyperspace, for r = 10 and k = 15, 20, 35, and 50 is presented as Box plots (Figure 3 4 A to D ). The OT had the highest average ED for any given k and r combination, followed by MOT. SU1000 and SU300 had the third and fourth highest distances. MM and SU1 were found to have, in general, the lowest ED s. In all cases the average ED s were of the same o rder of magnitude. OT, SU1000, and SU300 generally showed smaller interquartile ranges compared to MOT, MM and SU1. Results for r = 6, not presented, were similar to r = 10. All statistical tests for comparison of median ED resulted in rejection of the nul l hypothesis of equal medians at the 5% significance level. Parameter Screening Efficiency The average skill scores for all experimental combinations are presented in Table 3 4. For any given function configuration and r all sampling strategies had similar skill score values except for MOT. MOT had considerably lower skill scores for all G function configurations compared to other sampling strategies. On the other hand MOT performed relatively better than other sampling strategies for the G* functions. This indicated variation in the performance of MOT based on characteristics of the function or model being evaluated, potentially limiting its PAGE 91 91 utility across a range of models. Irrespective of the function, r , or sampling strategy the avg_ g scores were higher than avg_ R . Both the scores improved slightly as r increased from 6 to 10. Following Approach 1, SU300 had the highest AVG_ g and AVG_ R scores (Fig 3 5 A , B ). SU1000 had the second highest AVG_ g score , only slightly less than AVG_ g score for SU300 . On the o ther hand, MOT performed the worst with the lowest AVG_ g and AVG_ R scores. Overall, SU configurations performed better than the benchmark sampling strategies (OT and MOT) used in this study. The standard errors (Fig 3 5 A , B ) for all the sampling strategies were similar except MOT which had larger standard error compared to others. Results of parameter screening efficiency based on Approach 2 are presented in Figure 3 5 C and D . It can be seen that the sampling strategy with the lowest AVG R ank (the poorest performance) for both g and R scores was OT. While no particular sampling strategy was consistently ranked the best (lowest rank score) overall the SU strategies performed better than others with SU1000 and SU1 appearing in the first three best strategies for both scores and SU300 being the best strategy in the cases of R scores. Note that standard error for Approach 2 skill scores varied considerably from strategy to strategy with MOT showing the largest error for both scores indicating larger performance variability across function types. Discussion Theoretical Basis : Sampling for Uniformity Choosing the correct parameter probability distributions is critical in exhaustive uncertainty and sensitivity analyses ( Cacuci and Ionescu Bujor, 2004; Wallach et al ., 2006; Benke, 2008; Shin et al., 2013 ). Analogously, it is equally important for a sampling strategy to honor the probability distribution for which the sampling is being done at the minimum computational cost. A sampling strategy failing to generate the assigned distribution induces an extra layer of uncertainty which may lead to spurious results. Such a strategy should be avoided PAGE 92 92 irrespective of the type of numerical experiment, whether it be a full sensitivity analysis or a screening exercise. The re semblance between the theoretical and generated distributions for low dimensional problems (four parameters), the test cases used in Campolongo et al. (2007) and Ruano et al., (2012), is likely an artifact of the trajectory optimization i.e. distance maxim ization process. For example, for a single parameter model ED between trajectories will be maximum when levels are sampled evenly. These authors used visual inspection to draw conclusions about the supremacy of OT or MOT over MM pertaining to this criterio n. In theory, it is possible to generate samples such that distributions based on all points are perfectly discrete uniform. However, the number of trajectories required in such a case is likely to be high. Hence, as proposed by SU, it is logical to consid er only the first and the last points of the trajectories as each parameter takes exactly two values in a trajectory with the first and the last points having levels with constant offset . Statistical tests, such as the Chi Square test used in this study are usually applied at 0.05, though standard practice could vary from field to field. Here the generated parameter distributions were tested at three values of from 0.05 to 0.2. In spite of such loose criteria the average fraction of parameters fail ing the goodness of fit test were considerably large (10 70%) for MM, OT, and MOT with a general trend of higher number of failures with increasing number of parameters. From this perspective, SU performed the best with no failures which makes it a very in tuitive sampling strategy for models with a large number of parameter. Sample Generation Runtime Efficiency The main motivation behind the development of MOT was the impractically high sample generation time requirement of OT (Ruano et al., 2012) for mediu m to large problems. However, from the results here it was clear that there is not much advantage in implementing MOT if the PAGE 93 93 time requirement is the only selection criteria. For SU the time requirements can vary considerably depending on Q and r . With Q = 1 SU was the most time efficient strategy. Even with Q = 300, SU required about 20 minutes of computer time for k = 100, which can still be considered to be efficient and practically implementable. On the other hand, as seen in section 3.2, the sample gene ration time requirement for SU increased linearly with r as well as Q . For the r values used in this study SU configurations were time efficient . However , time requirements for SU will be equivalent or as expensive as OT and MOT if a large number of trajectories are generated ( r > 30) using a large number of repetitions ( Q > 500). Notice that the recommended number of trajectories reported in the EE literature vary from as few as 2 (Campolongo et al., 2011), to 10 50 (Campolongo et al., 2007), to r > 50 (Ruano et al., 2012), while a number of applications report r ~ 10 to be sufficient for parameter screening ( Francos et al., 2003; Cariboni et al., 2007; Carpena et al., 2007; Confalonieri et al., 2010; Kisekka et al., 2013 ). Herman et al. (2013b ) concluded that even with ~ 20 trajectories the relative importance of parameters can be identified well. Based on the above literature and experience of the authors, as few as ~10 trajectories should be sufficient for the screening experiments. Spread of Trajectories in Hyperspace and Parameter Screening Performance Given the intention of OT and MOT to cover a greater part of the parameter space by maximizing ED , one would expect these sampling strategies to be superior. However, this study indicated that the relationship between ED and parameter screening is not as strong as initially assumed in OT and MOT but has some significance. For example if the distance measure was the only criteria for better screening, OT should have performed better than MOT whi ch in turn should have performed better than SU300. On the other hand if spread was completely irrelevant then SU1 would have performed consistently better than SU300. However, these were not the PAGE 94 94 cases. This underlines the importance of the combined unifor mity and spread criteria. Identifying the exact role of trajectory spread in parameter screening should be a topic of further study. It is likely that a given sampling strategy is suitable for a particular type or types of models. As noted earlier, MOT was consistently the best sampling strategy for the G* configurations (Table 3 4) and the worst for the G configurations. Such knowledge about the model type sampling strategy relation may be useful for models with known behavior. However, prior notions abo ut the nature of newly developed models are rare and predicting the behavior of an existing model applied in a new set of conditions is difficult (Shin et al., 2013). Hence developing a sampling strategy that can effectively screen parameters for a wide ra nge of models is of paramount importance. There are also trade offs between the merits of a given sampling strategy. For example, OT produces the most widely spread parameter trajectories but takes the longest to generate the sample, while SU1 which needs the lowest generation time is not the best strategy in terms of screening. On the other hand, MM is considerably time efficient but performs poorly in producing uniform samples. In other words, no sampling strategy is the best or the worst . Instead, we desire a sampling strategy which performs gen erally well across all of the performance criteria. Table 3 5 presents a summary where each strategy has been assigned a qualitative ranking representing poor, acceptable and good performances. Considering approaches 1 and 2 separately while summing up the rankings enabled to assign double weightage to screening efficiency; the most important evaluation criteria. This qualitative comparison indicates that, generally, SU performs better than the benchmark trajectory based sampling strategies. PAGE 95 95 Other Consider ations In variance based GSA first order effects as well as total effects are calculated and studied. Campolongo et al. (2011) discussed the unification of screening and quantitative s such as Herman et al. unison. Based on this an argument, there is a need to maintain consistency across methods (screening and quantitative) in the consideration of sensitivity measures. For EE experiments it is a common practice to screen parameters based on * alone even though values are calculated. A study is required to evaluate the efficacy of as the surrogate for an interaction effect sensitivity index and effect of sampling strategy on this relation. The avg_ g scores for the OT strategy obtained by Campolongo et al. (2011) with the same B20, M20 and G*20 (referred as B, M, and G*10 , respectively, in their article) are similar to avg_ g scores in the case of the B function, somewhat different in the case of the G* function and considerably different in the case of the M function. One likely reason for these differences is the use of different parameter distribution truncation. As noted e arlier, we truncated parameter distributions at 12.5% and 87.5% following SimLab v2.2 for all sampling strategies, while Campolongo et al. (2011) used entire ranges. The use of entire parameter ranges may not prove to be an issue for test functions as all parameters had uniform distributions, but is not feasible for real world models with parameters distributions with infinite tails. Until now little attention has been paid to the effect of truncation on screening performance (Saltelli et al., 2009). Howeve r, with each unique model the corresponding suitable truncation could vary as well. Summary For high dimensional environmental models, parameter screening becomes an essential step before applying rigorous variance based GSA/UA techniques. The use of appr opriate PAGE 96 96 methodologies is of utmost importance when environmental models are used to support decision making and the selection of a suitable parameter sampling strategy is a multi faceted problem. In this paper we propose a new multi criteria trajectory bas ed parameter sampling strategy (SU) for calculating sensitivity measures based on the method of EEs for the parameter screening exercises. The main concept around which the SU is built is honoring the theoretical parameter distribution (i.e. discrete unifo rm distribution) along with good spread of parameter samples. While SU, in general, outperforms some of the trajectory based benchmark sampling strategies in most of the evaluation criteria (uniformity, spread, time requirements and screening efficiency) i ts implementation is flexible to control trajectory spread by selecting the desired value of the oversampling size. The future efforts to improve SU shall focus on sampling from a finer parameter grid ( p = 6, 8) and determination of optimal oversampling si ze. The advantage of SU as a sampling strategy is founded on combined sampling criteria. Exploration of other combinations of criteria e.g. such as other distance measures, will be helpful to the advancement of EE based parameter screening. Also, the stu dy of the relationship between EE based and GSA based interaction effect indices will be complementary to efforts focused on improvement, understanding, and potential unification of these sensitivity analysis approaches. PAGE 97 97 Table 3 1. Elementary effect metho ds and sampling strategies for parameter screening sensitivity analysis Method and Authors Availability Sensitivity Measures 1 Improvements Method of elementary effects or Morris Method (MM) (Morris, 1991) (1) As a part of SimLab v2.2 (2) Matlab ( http://ipsc.jrc.ec.europa.eu/?id=756 ) i and i Sampling strategy 2 , sensitivity measures Method of Optimized Trajectories (OT) (Campolongo et al., 2007) Matlab ( http://ipsc.jrc.ec.europa.eu/?id=756 ) i * and i Sampling strategy, sensitivity measures Combined LHS and one at a time sampling (van Griensven et al., 2006) NA 3 i * 4 Sampling strategy Simplex based sampling (Pujol, 2009) R ( http://cran.r project.org/web/packages/sensitivity/ ) i * and i Sampling strategy Cell based trajectory sampling (Saltelli et al., 2009) NA i * and i Sampling strategy Radial sampling (Campolongo et al., 2011) NA i * and i Sampling strategy Modified Optimized Trajectory Method (MOT) (Ruano et al., 2012) Matlab (personal communication with authors) i * and i Sampling strategy Sampling for Uniformity (SU) (This article) Matlab ( http://abe.ufl.edu/carpena/software/SUM orris.shtml ) i * and i Sampling Strategy 1 i = mean of distribution of elementary effects (EE) for i th parameter , i * = mean of distribution of absolute values of elementary effects (EE) for i th parameter , i : standard deviation of distribution of elementary effects (EE) for i th parameter 2 Original development 3 Not freely available 4 van Griensven et al. (2006) use sens itivity measure similar to i * PAGE 98 98 Table 3 2. Test function configurations and the number of important parameters (based on analytical expression/ previously calculated values of sensitivity indices) Function k B G G* M 1 O 2 i 3 #T 4 i #T i #T i #T i #T 15 NA 5 ~ 89% 4 ~ 10% 12 NA ~ 71% 5 20 ~ 34% 7 ~ 80% 7 ~ 24% 4 ~ 40% 10 NA 35 NA ~ 53% 10 ~ 79% 10 NA NA 50 ~ 79% 13 ~ 31% 11 ~ 52% 6 NA NA 80 ~ 10% 14 ~ 44% 13 ~ 65% 8 NA NA 100 ~ 52% 7 ~ 11% 18 ~ 36% 11 NA NA 1 2 Saltelli et al. (2008) 3 i sum of first order sensitivity indices 4 #T number of important parameters based on total effect sensitivity indices ( ST i s). The parameter was considered important if ST i > 5% 5 Not Applicable PAGE 99 99 Table 3 3. The average fraction of parameter distributions generated by the six sampling strategies failing the Chi square test for discrete uniform distribution at three significanc e levels. Sampling strategies are: Morris Method (MM), Optimized Trajectories (OT), Modified Optimized Trajectories (MOT), Sampling for Uniformity (SU). MM OT MOT SU 1 2 k 3 0.2 0.1 0.05 0.2 0.1 0.05 0.2 0.1 0.05 0.2 0.1 0.05 15 0.30 0.73 0.73 0.08 0.57 0.57 0.19 0.66 0.66 0.00 0.00 0.00 20 0.34 0.75 0.75 0.11 0.59 0.59 0.22 0.67 0.67 0.00 0.00 0.00 35 0.34 0.75 0.75 0.16 0.66 0.66 0.25 0.70 0.70 0.00 0.00 0.00 50 0.35 0.76 0.76 0.20 0.65 0.65 0.26 0.72 0.72 0.00 0.00 0.00 80 0.35 0.76 0.76 0.24 0.69 0.69 0.28 0.72 0.72 0.00 0.00 0.00 100 0.28 0.72 0.72 0.25 0.69 0.69 0.28 0.72 0.72 0.00 0.00 0.00 1 Results for all Sampling for Uniformity (SU) configurations i.e. SU1000, SU300 and SU1 were identical. 2 is the significance level at which statistical tests were performed 3 k is the number of parameter used for sample generation PAGE 100 100 Table 3 4. Summary of screening efficiency skill scores (avg_g and avg_R ) for six sampling strategies: Morris Method (MM ), Optimized Trajectories (OT), Modified Optimized Trajectories (MOT) and Sampling for Uniformity in three configurations (SU1000, SU300, and SU1) Test Function avg_ g avg_ R r MM OT MOT SU1000 SU300 SU1 MM OT MOT SU1000 SU300 SU1 G15 6 0.9980 0.9980 0.9870 0.9970 0.9965 0.9955 0.8145 0.7800 0.1990 0.8130 0.8120 0.8040 10 0.9995 0.9990 1.0000 0.9995 0.9995 0.9985 0.9135 0.8880 0.3835 0.8975 0.9015 0.8800 G*15 6 0.8835 0.8828 0.8898 0.8773 0.8750 0.8772 0.0992 0.0857 0.2325 0.0970 0.0968 0.0953 10 0.8957 0.9002 0.8998 0.8858 0.8890 0.8887 0.0993 0.1123 0.2382 0.1010 0.1018 0.1005 O15 6 0.8614 0.8586 0.8769 0.8689 0.8666 0.8609 0.1954 0.2111 0.2049 0.2171 0.2174 0.2151 10 0.8754 0.8897 0.8963 0.8946 0.8940 0.8983 0.2143 0.2314 0.2320 0.2406 0.2426 0.2369 B20 6 0.9810 0.9920 0.9929 0.9824 0.9848 0.9826 0.2537 0.3124 0.3549 0.3010 0.2901 0.2517 10 0.9954 0.9993 0.9986 0.9977 0.9940 0.9957 0.2949 0.3596 0.3759 0.3199 0.3337 0.2889 G20 6 0.9746 0.9710 0.8767 0.9720 0.9704 0.9726 0.3864 0.3820 0.0824 0.3963 0.3849 0.3963 10 0.9901 0.9857 0.9191 0.9859 0.9877 0.9877 0.4981 0.4624 0.1010 0.4693 0.4926 0.4611 G*20 6 0.8415 0.8420 0.8383 0.8498 0.8405 0.8245 0.2450 0.2688 0.3243 0.2303 0.2510 0.2300 10 0.8783 0.8775 0.8888 0.8805 0.8883 0.8630 0.2725 0.2713 0.3748 0.2370 0.2450 0.2610 M20 6 0.8696 0.8546 0.8685 0.8587 0.8695 0.8720 0.1524 0.1511 0.1208 0.1656 0.1642 0.1550 10 0.9090 0.8998 0.9100 0.9135 0.9073 0.9154 0.1837 0.1721 0.1398 0.1786 0.1858 0.1874 G35 6 0.9957 0.9958 0.8900 0.9942 0.9968 0.9960 0.1796 0.1690 0.0234 0.1712 0.1746 0.1822 10 0.9991 0.9996 0.9290 1.0000 1.0000 0.9972 0.1986 0.2044 0.0200 0.1614 0.1890 0.1914 G*35 6 0.9974 0.9986 1.0000 1.0000 0.9974 0.9990 0.2489 0.2648 0.4720 0.2544 0.2692 0.2724 10 1.0000 0.9998 1.0000 1.0000 1.0000 1.0000 0.2972 0.3026 0.5082 0.3542 0.3222 0.3196 B50 6 0.7857 0.7908 0.8077 0.8015 0.7943 0.7961 0.0685 0.0615 0.0846 0.0672 0.0663 0.0682 10 0.8077 0.8277 0.8231 0.8222 0.8218 0.8122 0.0768 0.0585 0.0769 0.0714 0.0637 0.0743 G50 6 0.8991 0.8727 0.6727 0.9016 0.8984 0.8931 0.1193 0.1164 0.0318 0.1000 0.0951 0.1154 10 0.9256 0.9236 0.7182 0.9231 0.9269 0.9218 0.1235 0.1382 0.0227 0.1207 0.1124 0.1277 G*50 6 0.6323 0.6067 0.8500 0.6289 0.5947 0.6107 0.2363 0.2000 0.4333 0.1789 0.1593 0.1945 10 0.6843 0.6467 0.8417 0.6550 0.6393 0.6612 0.2667 0.1667 0.4500 0.2257 0.2160 0.2457 B80 6 0.7379 0.7266 0.7387 0.7380 0.7399 0.7317 0.1194 0.1089 0.1046 0.1151 0.1136 0.1059 10 0.7827 0.7816 0.7804 0.7951 0.8073 0.7714 0.1384 0.1390 0.1363 0.1250 0.1263 0.1299 G80 6 0.8522 0.8495 0.5191 0.8265 0.8189 0.8383 0.1069 0.1046 0.0246 0.0808 0.0882 0.1069 10 0.8832 0.8832 0.5914 0.8832 0.8802 0.8683 0.1191 0.1298 0.0294 0.0717 0.0677 0.1042 G*80 6 0.8553 0.8740 0.9228 0.8398 0.8628 0.8515 0.2075 0.2005 0.4368 0.1068 0.1755 0.1775 10 0.8805 0.8868 0.9390 0.8813 0.8765 0.8660 0.2525 0.2563 0.4940 0.2478 0.2913 0.2410 PAGE 101 101 Table 3 4. Continued Test Function avg_ g avg_ R r MM OT MOT SU1000 SU300 SU1 MM OT MOT SU1000 SU300 SU1 B100 6 0.9423 0.9429 0.9040 0.9729 0.9706 0.9517 0.2100 0.2009 0.2437 0.1629 0.2120 0.2143 10 0.9806 0.9803 0.9520 0.9786 0.9929 0.9814 0.2854 0.2569 0.2603 0.2129 0.2543 0.2703 G100 6 0.9152 0.8999 0.7048 0.8899 0.9187 0.9001 0.0642 0.0629 0.0256 0.0800 0.0738 0.0477 10 0.9329 0.9334 0.7844 0.9311 0.9238 0.9232 0.0673 0.0597 0.0166 0.0622 0.0816 0.0584 G*100 6 0.8787 0.8813 0.9576 0.8415 0.8545 0.8453 0.1245 0.0956 0.3331 0.0942 0.0845 0.0958 10 0.9109 0.9238 0.9673 0.9091 0.9091 0.8962 0.1282 0.1444 0.3564 0.0909 0.1115 0.1280 PAGE 102 102 Table 3 5. Qualitative assessment 1 of sampling strategies used in this study across all evaluation criteria Sampling Strategy Evaluation Criteria Total Score Uniformity Time Spread Screening Efficiency Approach 1 Approach 2 MM * *** * *** * 9 OT ** * *** * * 8 MOT * * *** * ** 8 SU1000 *** ** ** ** *** 12 SU300 *** ** ** *** *** 13 SU1 *** *** * ** *** 12 1 * poor, ** acceptable, *** good PAGE 103 103 Figure 3 1. Schematic representation of parameter levels and constant offset coordinates/ jump size for an arbitrary parameter in A ) unit parameter space and B ) truncated unit parameter hyperspace A B PAGE 104 104 Figure 3 2. Schematic representation of sample generation using Sampling for Uniformity (SU) for the case of 2 parameters, 4 trajectories A ) Step 1, B ) Step 2, and C ) Step 3. F, I and L correspond to the first, intermediate and the last point of each trajectory, respectively. Numbers following the letters F, I and L indicate the trajectory numbers. A B C PAGE 105 105 Figure 3 3 . Time required for sample generation of size k at two trajectory configurations ( r = 6 and r = 10) using six parameter sampling strategies: Morris Method (MM) with Q = r (by definition), Optimized Trajectories (OT) with Q = 1000, Modified Optimized Trajectories (MOT) with Q = 1000, Sampling for Uniformity with Q = 1000 (SU1000), Sampling for Uniformity with Q = 300 (SU300), and Sampling for Uniformity with Q = 1 (SU1) PAGE 106 106 Figure 3 4. Box plots of the Euclidean distances (measure of trajectory spread): Morris Method (MM), Optimized Trajectories (OT), Modified Optimized Trajectories (MOT), Sampling for Uniformity with Q = 1000 (SU1000), Sampling for uniformity with Q = 300 (SU300) and Sampling for Uniformity with Q = 1 (SU1); across range of model sizes (i.e. number of parameters, k ) for r = 10 and for A ) k = 15, B ) k = 20, C ) k = 35, D ) k = 50. A B C D PAGE 107 107 F igure 3 5. Parameter screening efficiencies for six sampling strategies: Morris Method (MM), Optimized Trajectories (OT), Modified Optimized Trajectories (MOT), Sampling for Uniformity with Q = 1000 (SU1000), Sampling for Uniformity with Q = 300 (SU300), a nd Sampling for Uniformity with Q = 1 (SU1); across all test functions using Approach 1 ( A , B ) and Approach 2 ( C , D ) A B C D PAGE 108 108 CHAPTER 4 GLOBAL SENSITIVITY AND UNCERTAINTY ANALYSIS OF THE WATERSHED ASSESSMENT MODEL (WAM) Sensitivity analysis (SA) is an essential ingredient of the modeling process (Saltelli et al., 2004) and becomes critical when the model is very compl ex with hundreds of parameters, a typical attribute of large scale hydrologic and water quality (H/WQ) models (Gan et al., 2014) . The outcomes of a sensitivity analysis help the calibration validation process. The need for Carpena et al., 2006; Shirmohammadi et a l., 2006). From the literature survey presented in Chapter 1 it is clear that a considerable amou n t of attention is being paid to sensitivity and uncertainty analysis (UA) studies in H/WQ modeling with a large number of articles on th ese topic s published i n top journals ; such as Environmental Modelling an d Software, Hydrological Processes, Journal of Hydrology, Water Resource Research etc., during the last 10 15 years. At the same time it i s evident that SA/UA of H/WQ models is still an evolving field. Earl ier applications of SA/UA to large scale H/WQ models were restricted by the availability of computational technology. The r ecent advent in para llel computing has enabled mode lers to apply the most computationally expensive and the most robust SA/UA techniques. In spite of this , the practice s of using local SA techniques (see Table 1 5) as well as reducing model dimensionality based on previous model analysis in different conditions have continued. Both these practi ces are scientifically inappropriate from the large scale H/WQ modeling perspective as (a) these models are so complex that they are likely to involve non linearities and parameter interactions which cannot be detected by local SA and (b) any new applicati on of a model is unique and it is difficult to predict which model parameters are influential or non influential based on prior model analysis in other environment s . PAGE 109 109 The Watershed Assessment Model (WAM) is a watershed scale simulation model that rep resents terrestrial hydrology (Bottcher et al., 2012). The s alient features of WAM were presented in Table 1 1. An e arlier version of WAM , i.e. WAMView , was one of the models recommende d by the US EPA for Total Maximum Daily Load (TMDLs) and watershed protection planning ( US EPA, 2013) ( http://www.epa.gov/athens/wwqtsc/index.html ) . WAM has been successfully implemented in a number of watershed planning and management projects throughout the state of Florida during the last two decades (SWET, 2014) , a few of which are as listed below. Lake Okeechobee Watersheds for P budge analysis ( various projects since early 2000 ) ( Zhang et al., 2006; Keener et al., 2007; SWET, 2009a; SFWMD, 2010; Chebud et al., 2011 ). Modeling the Upper Peace River and three lakes for TMDL development ( for FDEP ) (SWET, 2009 b ) St. Johns River sub watersheds for the assessment of land use scenarios on hydrology and water quality ( late 1990s ) (SWET, 1997 ) Suwan ee River watershe d assessment for water quality (mid 1990s) and BMP evaluations ( 2008 ) (SWET, 201 4 ) TMDL development for the Mayakka and Hillsborough rivers ( for US EPA ) FDEP (SWET, 201 4 ) Conceptual design and testing of regional stormwater treatment areas and conditions modeling (SWET, 201 4 ) It can be argued that it is one of the most suitable model for Florida due to (a) WAM uses field scal e models that can handle divers e hydrologic conditions that are typical to the state of Florida and ( b ) an extensive model database for Floridian conditions (soils, BMPs etc.) that has been developed over the period of two decades from 1995 201 4 . Note that among the most popular watershed models only a few (HSPF and MIKE SHE) have been successfully applied in PAGE 110 110 Florida (Bergman et al., 2002; Murray et al., 2004; Ross et al., 2005; Hosseinipour, 2006; Jaber and Shukla, 2007; Kang and Gilbert , 20 13 ). In spite of its suitability as a watershed model in general WAM has several weaknesses e of the major issues identified by the independent peer review committee for WAM (Graham et al., 2009) include d ( a ) lack of transparency in model documentation about model assumptions, limitations , and methodologies , ( b ) no formal sensitivity analysis in any model application , ( c ) lack of thorough model evaluation criteria i.e. model calibration validation procedure , and ( d ) need for the predictive uncertainty in model applications. As a result, WAM underwent a rigorous model documentation process in the l ast four years (2010 2014) to address the recommendations provided by the peer review committee. The p roced ure for WAM calibration and validation has been published in Bottcher et al. (2012). However, issues related to model sensitivity and uncertainty hav e not been studied t o date. Such studies are essential to better understand the WAM processes, its limitations, and the overall reliability of the model. The Numeric Nutrient Criter i a for Florida (US EPA, 2010) makes evaluation and implementation of various water resources policies (e.g. TMDLs) imperative. Due to its application history in the state , WAM is likely to be the first choice model for new watershed management studies . This further increase s the need to focus on a critical aspect reliability studies of WAM . In t his work on SA of WAM we have focused on B asin Land Area to Stream Routing ( BLASROUTE ) i.e. the flow and water quality constituent routing sub module of WAM . The re asons for considering only the routing part include: ( a ) a number of Basin Unique Cell Shell ( BUCHSHELL ) (source scale simulation sub module) parameters are hard coded i.e. not available for user to adjust, ( b ) simplified approaches and empiricism a re more pronounced in PAGE 111 111 BLASROUTE (Graham et al., 2009) , ( c ) almost all of the BLASROUTE parameters can be and need to be adjusted by the user for calibration purpose s , ( d ) even for a small and simple watershed (watershed with very few types of soils and land uses) the total adjustable parameters (BUCSHELL and BLASROUTE parameters together) that need to be considered in SA/UA will be of the order of 10 3 10 4 making these a nalyses extremely time and resource consuming. The specific objective s of this part of this research are to perform (a) global SA of WAM (b) UA of selected WAM o utputs, and (c) parameter estimation through multi objective Monte Carlo Filtering (MCF). For the EE analysis sampling strategy developed in Chapter 3 is applied . The Taylor Creek Nubin Slough i.e. S 191 watershed has been used as a test case for WAM appl ication. Material and Methods Model WAM is a GIS based H/WQ model that represents water quantity and quality responses within the terrestrial portion of the hydrologic cycle based on detailed characterization data (Bottcher et al., 2012). WAM simulates th e constituents important to eutrophication processes in water bodies (water, total suspended solids, biological oxygen demand, soluble and particulate forms of nitrogen and phosphorus) within a watershed. WAM shares some of its features with other popular watershed models indicating its soundness as a watershed model in general (Table 1 1 ) . It has two main sub modules BUCSHELL and BLASROUTE which generate flow/chemical constituents and route them through the watershed/stream network, respectively. These su b modules are run consecutively and independently one after the other. PAGE 112 112 Source cell simulations Before the physical processes are simulated, the watershed is gridded into a rectangular array based on a user defined cell size (typically 1 ha) as shown in Fig ure 4 1 . Cells that have identical characteristics are termed a unique cell type and are simulated only once to reduce the simulation time . The BUCSHELL program simulates the flows and water quality constituents from each cell on a daily basis using one of the three field scale models (GLEAMS, EAAMOD or Special Case). BUCSHELL selects field scale models based on land use type and soil characte ristics (Figure 4 2), allowing for a broad spectrum of conditions to be simulated. GLEAMS ( Knisel, 1993 ) and EAAMOD ( B ottcher et al., 1998 ; Zhang and Gornak, 1999 ) are physically based field scale models for agricultural land uses and have a number of par ameters. The third field scale model, termed as the Special Case ( SWET, 2008; Bottcher et al., 2012 ) , is designed to simulate land uses that cannot be represented by typical crop models, such as mines, open water, and wetlands. It uses a simplified water b alance method to calculate daily water volumes contributing to surface runoff and groundwater. A user defined nutrient concentration is applied to runoff and groundwater percolation when it occurs. Flow and nutrient routing Water and nutrient s generated on source cells by the BUCSHELL are routed to their ultimate destinations using the BLASROUTE sub module. The paths and distances flows and flow constituents need to travel are determined by source cell characteristics , topography, and the features (wetlands /depressions) they come across enroute to corresponding model stream reaches. WAM uses six different distance grids for source cell to stream routing as shown in Figure 4 3. Constituents are routed using overland flow and groundwater flow unit hydrographs , delay factors, and generalized attenuation model s . Once flows and loads have been delivered to the stream reaches, they are routed through the network to the basin outfall using a modified PAGE 113 113 linear reservoir routing technique developed by SWET ( Jacobson and Bottcher, 1998 ) for solving the equations of unifor . This in stream routing technique can handle hydraulic loops and structures. In stream routing and attenuation are affected by channel slope, cross section profile, and features such as gates, weirs etc. within the stream network. A ttenuation models for overland flow, groundwater flow, streamflow , and cell to reach flow routing are represented by eq s . 4 1 , 4 2, 4 3, and 4 4 , respectively . (4 1) (4 2) (4 3) (4 4) where C = outflow concentration (ppm) , C o = inflow concentration (ppm) , C b = background concentration (ppm) , q = discharge (m 3 s 1 ) , d = flow distance, a and b = empirical parameters , = time step (s), R = hydraulic radius (m) , T = time of travel to the feature (d), v = flow velocity (m d 1 ), k = constant time delay coefficient (d), U_HYDRO = unit hydrograph . Study Watershed and Model Set up For this study we used Taylor Creek Nubin Slough watershed (Hydrologic Unit Code 03090102) in Okeechobee County , Martin County, and St . L ucie County, Flor ida. This watershed, located on the north east shore of Lake Okeechobee, drains ~ 485 km 2 of land in to the lake through the water control structure (spillway) S 191 . The t opography is generally flat and the main direction of the slope is from the eastern b oundary to the western boundary of the watershed (Figure 4 4 A ). Taylor Creek drains the northern part of the watershed . Nubin Slough, Mosquito Creek, Henry Creek, and Lettuce Creek drain the central and the southern part s of the PAGE 114 114 watershed to the l e vee L 63 which is ultimately drained by S 191. C hannels are generally poorly defined, have broad flood plains, and the streamflow is slow moving (Keener et al., 2007). Dairy and beef cattle farming is the ma jor activity in the watershed. Animal grazing pasture s (improved ~ 52% , intensive ~ 4%, unimproved and other types ~ 5% ) are located throughout the watershed accounting for ~61% of the total area (Table 4 1, Figure 4 4 B ) . Among other land uses wetlands (~9%) and agriculture (other than pastures) (~9%) form substantial portions of land cover . The dominant soil classes in the watershed are Immokalee and Myakka covering ~27% area (each) (Figure 4 4 C ) . Overall the soils are poorly drained, coarse textured soils with low phosphorus retention capacity. This water shed has been identified as one of the main contributors of phosphorus load to the lake ( McCormick et al., 2010; Zhang and Sharfstein, 2013 ). This has led to a nu mber of projects such as Sormwater Treatment Ares (STAs) , a gricultural best management practic es (BMPs), wetland restoration programs, etc. ( SWET, 2009 a ; S F W MD, 2010 ). Over the years WAM has been implemented in a number of modeling studies for the Lake Okeechobee watersheds including the S 191 basin ( e.g. Keener et al., 2007; SWET, 200 9a; SFWMD, 2010 etc. ). The model used as a test bed for this work is the same as the base scenario used in the WAM Documentation and Validation Project ( SWET, 2011 ). The p reviously calibrated model set up was obtained from SWET Inc. (personal communication). BUCSHELL was simulated once to generate daily time series of flows and water quality constituent loads from individual model cells. These day files were used as deterministic model input to BLASROUTE in this study. The m odel was simulated for three years from 20 02 2004 which was the calibration period in the WAM Documentation and Validation project (SWET, 2011) . The BLASROUTE parameter s can be categorized into six types (a) water quality constituent PAGE 115 115 attenuation (31 for each wetland type and dryland , 18 for each r each type ), (b) unit hydrographs (runoff and subsurface), (c) delay factors ( 4 for each dryland, wetland, and depression), (d) channel cross section (c/s) geometry (c/s profile for each calculation reach) , (e) channel hydraulic properties and leakage (for each calculation reach) , (f) boundary characteristics (tidal boundaries), and (g) hydraulic structure related parameters. Since S 191 ha d 2 types of reaches, 7 types of wetlands, no flow control structures (except the outlet), and no tidal boundary t he total number of parameters w as 297. Parameter details are presented in Appendix C. Evaluation Methodologies In this study WAM/BLASROUTE has been evaluated using three techniques namely, (a) global sensitivity analysis, (b) uncertainty analysis, and (c) parameter estimation. Details of these methodologies are presented in following sections. Global s ensitivity analysis The SA of WAM/BLASROUTE was performed by sequential application of th e method of EE (Morris, 1991) Saltelli (2002) . This two step SA framework has been recommended in a few H/WQ modeling S A studies ( e.g. oz Carpena et al., 2007 ) . The EE analysis i s a low cost qualitative method and is used to identify influential model parameters which are then consider ed for the qu a Some of the recent studies e.g. Foglia et al. (2009), Srivastava (2013), Srivastava et al. (2014) have suggested that considering the impractically large computational resources and time requirements of robust SA techniques (variance based SA) . However, variance based SA techniques prov ide more insights about the model and should be used whenever possible. The f ollowing sections provide details of the SA techniques used in this study. PAGE 116 116 The method of EE : The first step in the EE method is to generate a parameter sample which is done from t he k dimensional unit hyperspace divided into a p level grid. A p level grid implies that each factor can take only p distinct values {0, 1/( p 1), 2/( p r trajectories (or paths), each consisting of k +1 points are generated . In other words , the total number of model runs required for the EE analysis is r ( k +1) . Within any given trajectory two samples differ from each other in p /[2( p 1)]. The EE associated with the i th model parameter is calculated as: ( 4 5 ) where y is the model output to the corresponding factor set, x i is the i th input factor. r elementary effects are obtained from the entire sample for each of the parameter s . The d istribution of EE for the i th parameter is analyzed statistically to calculate sensitivity measures the mean and standard deviation of the EE distributions ( i and i ) . Sensitivity measure i suggested by Morris (1991) w as found to be inefficient and w as subsequently modified by Campolongo e t al. (2007). Th is new measure , the mean of absolute values of EEs ( i * ) ha s been found to produce results consistent with the total effect sensitivity index in variance based quantitative SA . For each model output of interest parameters are plotted in i * i space. The h igher the value of * of a parameter the higher is its overall influ ence on that output. Parameters that are well separated from the origin of * space are important while the ones close to the origin are less important and can be omitted from more detailed analyses . In this study we have calculated both and * and used and * spaces to plot EE analysis results in order to (a) get more i nsights of EE analysis and (b) identify flaws in WAM coding base d on sign of (discussed later) . PAGE 117 117 The method of : , 1993) is a variance based global sensitivity analysis technique. Two sensitivity indices are calculated per parameter. The first order sensitivity index ( S i ) of a parameter indicates the contribution of that parameter to total output variance ( V ( Y )) while the total effect sensitivity index ( ST i ) is its overall contribution (i.e. contri bution including its interactions) to total output variance. ( 4 6 ) ( 4 7 ) where is the variance of expectation of Y over the i th parameter , is the variance of expectation of Y over all but the i th paramete r and similarly is the expectation of variance of Y over all but the i th parame ter . Other d etails about this method are presented in Chapter 1 where it was also noted that the required N (2 k +2) model simulations w ith N > 500 is a recommended value. Such a large requirement of model runs was due to the use of two resampling matrices . Saltelli (2002) showed that it is possible to calcu late both indices at a much reduced cost of N ( k +2) with the use of only one resampling matrix. In Saltelli et al. (2010) it was found that among several estim ates of variance decomposition the one proposed by Jansen (1999) works the best. Following these developments the Saltelli (2002) scheme for sampling and Jansen (1999) formula for variance decomposition. Uncertainty analysis and parameter e stimation Uncertainty analysis of WAM/BLASROUTE outputs of interest was performed by plotting probability density function s (PDF s ) , cumulative density function s (CDF s ) , and PAGE 118 118 calculating the 95% confide nce interval (2.5% 97.5%) . Model runs f rom were used for this purpose. Parameter estimation was performed by MCF which categorizes model runs into two classes based on user defined criteria. This technique is essentiall y the same as the Generalized Likeli hood Uncerta inty Estimation (GLUE) method o logy (Beven and Binley , 1992) if the output of interest is a goodness of fit/likelihood measure. In this study a multi criteria MCF was performed. Details of model outputs and criteria use d are described in the experiment deta ils section. Experiment Details Parameter distributions Parameter ranges and their distributions can have substantial effect on SA results ( Cacuci and Ionescu Bujor, 2004; Wallach et al., 2006; Benke, 2008 ). Various methods to determine parameter distribut ions include (a) actual measurements, (b) literature review, and (c) expert opinion etc . However, it is not always possible to obtain detailed distribution s . Haan et al. (1995), oz Carpena et al. (2007) etc. have suggested the use of simple distributions such as uniform (when only the parameter range is known) or triangular (when a typical value is known along with the range). Other distributions can be assigned if additional information is available. If only a typical value is know n then the parameter can be assigned a uniform distribution with +/ 20 25% variation ( Zajac, 2010 ; Wallach et al., 2014 ). Among the BLASROUTE water quality attenuations parameters , a b are empirical in nature while C b an be estimated from monitoring data. a b parameters were assigned l og10 u niform (logarithm of values to base 10 being distributed uniformly) and u niform distributions based on the expert opinion of WAM developers and default parameter values i n the WAM manual (SWET, 2014) C b PAGE 119 119 reach types were estimated from the measured data at USGS gages in S 191 (USGS , 2014 ) ( http://waterdata.usgs.gov/nwis/dv/ ), literature review ( Byrne and Wood, 2012 ) and were assigned u niform distributions. On the other hand C b parameters for wetlands and dryland were assigned u niform and l og10 u niform distributions based on limited literature ( Kadlec and Knight, 1996; Wallace and Knight, 2006; Kadlec and Knight, 2009 ). For runoff/surface and subsurface hydrographs 4 profiles were created for each type and were assigned discrete uniform distribution . These profiles are shown in Figure C 1 and C 2 (Appendix C) . analysis 16 additional profiles of the runoff unit hydrograph were created (Figure C 3) . Time delay parameter distributions were estimated from discussion with the WAM developers. The (wid_err) (dep_err) are surrogate parameters which are not present in WAM. They were used as a proxy for the chann el geometry model in WAM . WAM developers suggest users to obtain reach profiles from actual measurements (SWET, 2014) and use them instead of the c/s geometry model . Since a previously calibrated set up for the S 191 was used instead of regenerating channel profiles for each simulation using the c/s profile model , it was decided to change the width and depth of each calculation reach indirectly by apply i ng percentage error parameters. It was assumed that these two c/s error parameters were distributed uniformly between 20% to +20%. roughness coefficient ( n ) was assigned a wide range of 0.02 to 0 .08 and a u niform distribution based on value s reported for natural channels and constructed canals with varying degrees of bed properties ( Chow, 1959; Florida Department of Transport, 2004 ). Details of all S 191 model parameters are presented in Table C 1 (Appendix C). O utputs of interest In most of the WAM applications to the Lake Okeechobee watersheds the modeling objective was to estimate annual flows and nutrient loads (SWET, 2009a; McCormick et al., PAGE 120 120 2010; SFWMD, 2010; Zhang and Sharfstein, 2013). Considering these past applications and the fact that WAM has not been tested for sensitivity before this study, a total of 11 model outputs for four water quality constituents ( T otal Phosphorus TP, T otal Nitrogen TN, Biochemical Oxygen Demand BOD, and Total Suspended Solids TSS) and discharge wer e analyzed using the method of EE to determine important BLASROUTE parameters using the Nash Sutcliffe Efficiency (Nash and Sutcliffe, 1970) as goodness of fit measure. Later , SA and UA it was decided to omit BOD and TSS related outputs. This was done to reduce the number of model runs required by SA/UA. Since BLASROUTE processes for water quality constituents were found to be independent (discussed later) omission of BOD and TSS parameters would not have affected TP and TN outputs. Additi onal l y , for the purpose of multi criteria MCF , an other goodness of fit measure /likelihood measure ; Percentage Bias (PBIAS) was also considered . M athematical expressions for NSE and PBIAS are as presented in eqs 4 8 and 4 9 . A s ummary of outputs considered for each of the analyses is presented in Table 4 2. Table 4 3 summarizes the criteria used in MCF which were based on the recommendation of Moriasi et al. (2007). (4 8 ) (4 9) where O i = i th observed value, i = average of n observations, P i = i th predicted/simulated value, n = number of observation prediction pairs. NSE varies between and 1. 1 indicates a perfect model while NSE < 0 indicates the me an of observations is better than the model. A PBIAS PAGE 121 1 21 value of a perfect model is 0. A n egative PBIAS indicates the model overe stimates while p o s itive value indicate s the model under e stimates observations. Other details For the EE analysis 10 parameter trajectories were generated (i.e. r = 10) using the Sampling for Uniformity (SU) sampling strategy developed in C hapter 3 of this work with 4 parameter levels and 300 repetitions . A pproximately 7 hours were required to generate this sample for 297 parameters . For comparison (and an additional evaluation of the SU strategy), sample generation using the method of OT (Campolongo et al., 2007) , with an oversampling size of 500 (lower end of the recommended oversampling size range) , took ~ 28 hours. Each model simulation required ~ 350 sec on a Intel Core 2 Duo (3 GHz) processor with 8 GB RAM desktop machine resulting in ~ 290 computational hours for the entire EE experiment i.e. 10x(297+1) = 2980 simulations. with 27 paramet ers identified in the EE analysis, N was set to 1000 which resulted in 1000 x (27+2) = 29000 model runs. The s analysis was generated using the R Sensitivity Package (Pujol, 2014). analysis were performed using the High Performance Computing Facility of the University of Florida (UF HPC) ( http://www.hpc.ufl.edu/ ). Results Global Sensitivity Analysis Elementary Effect sensitivity analysis Results of EE analysis for 11 WAM/BLASRO UTE outputs are shown in Figures 4 5 to 4 9 . Plots on the left side of these figures show EE results in the * space while the ones on right hand side show results in the space. Parameters for which the * value was not equal to the absolute value of corresponding we re indicated in blue. Also, on the plots two lines PAGE 122 122 corresponding to +/ 2 SEM ( standard error of the mean) we re plotted. Morris (1991) proposed that parameters located outside the wedge formed by the two lines can be regarded as important. For two of the three flow related outputs, NSE daily flow and NSE monthly flow, n (Fig 4 5 A to D ) . D_GW_vel (dryland to gr oundwater velocity) was found to be the most important factor for the t otal flow volume while it was the third most important parameter for NSE monthly flow. Depth error (dep_err) was found to be an important parameter for all the three outputs. EE analysis results for the TP re lated outputs NSE monthly avg . TP con . , NSE monthly TP load, and t otal TP load are shown in Figures 4 6 a to f. D_Sl_P_a and D_Sl_P_b i.e. dryland soluble P parameters a and b were found to be important parameters for all outputs . St_Sl_P_a ( stream soluble P parameter a ) was found to be important for load related outputs while St_Sl_P_Cb (stream soluble P background concentration) was found to be important for monthly average c oncentrations. n was considerably important for NSE monthly TP load. In this work TN was calculated as the sum of organic nitrogen, nitrate, and ammonia as WAM does not calculate nitri t es. Sensitivity analysis results for the TN outputs (Figs 4 7 A to F ) indicated that D_Sl_ON_a (dryland soluble organic nitrogen parameter a ) and D_Sl_ON_b (dryland soluble organic nitrogen parameter b ) are the first and the second most important parameters fo r all the three TN outputs. St_Sl_ON_a (stream soluble organic nitrogen parameter a ) had substantial impact on NSE monthly concentration but was not found to affect the other two outputs. For NSE monthly TN load D_GW_vel (dryland to groundwater velocity) , D_Sl_ON_Cb (dryland soluble organic nitrogen background concentration) , FM_ON_Cb PAGE 123 123 (freshwater marsh organic nitrogen background concentration) , D_Sd_ON_a (dryland sediment organic nitrogen parameter a ) were found to be important as well. For the t otal BOD load D_Sl_BOD_Cb (dryland soluble BOD background concentration) , D_Sl_BOD_a (dryland soluble BOD parameter a ) , D_Sl_BOD_b (dryland soluble BOD parameter b ), and FM_Sl_BOD_Cb (freshwater marsh soluble BOD background concentration) were the most important parameters (F ig 4 8 A , B ). Similarly for t otal TSS load D_TSS_Cb (dryland TSS background concentration) , D_TSS_a (dryland TSS parameter a ) , D_TSS_b (dryland TSS parameter b ) , and FM_TSS_Cb (freshwater marsh TSS background concentration) were identified as important (Fig 4 9 A , B ). F or a number of outputs clusters of parameters were observed somewhat separated from the origin of * plots (as well as ) but not as far away as important parameters were . Almost all t hese less importan t parameters (total 27 parameters) analysis considering the qualitativeness of the EE analysis as found in Chapter 3 and recommendation of WAM developers . Table 4 4 and 4 5 summarize these parameters. flow, NSE monthly flow, and t otal flow volume are shown in F igure 4 10 A , B , and C , respectively. For NSE daily flow D_RO_kDF explain ed 91% of the total output variance by itself ( i.e. S i = 91% ) while its overall contribution (including interactions) was 97% (i.e. ST i = 97%) . For all other parameters S i s as well as ST i s we re extremely small. In the case of NSE monthly flow D_RO_kDF, n, RO _HYD (runoff unit hydrograph), D_GW_vel, and dep_err contribute d (total effect) 56%, 29%, 22%, 13% , and 10% to the total variance , respectively. Their respective S i values we re 27%, 11%, 6%, 12%, and 6% indicating involvement of PAGE 124 124 D_RO_kDF, RO_HYD, and n in interactions. On the other hand for t otal flow volume D_GW_vel was the most dominating parameter with S i ~ 85% and ST i ~ 87%. Figure 4 11 A indicates St_Sl_P_a, D_Sl_P_a, D _Sl_P_ b , and St_Sl_P_Cb were the most importa nt parameters w ith ST i s ~ 57%, 43%, 24%, and 16%, respectively. Also, these parameters were involved in considerable interactions as the differences between corresponding S i and ST i were ~ 28%, 33%, 15%, and 12%, respectively. D_Sl_P_a was the most important parameter for NSE monthly TP load and t otal TP load with ST i ~ 71% & 60% and S i ~ 42% & 51% , respectively ( F igs 4 11 B and C ) . St_Sl_P_a and D_Sl_P_b are the other two important parameters for both NSE monthly TP load and t otal TP lo ad S i s and ST i s for the two outputs were 12% & 21% and 20% & 29%, respectively. corresponding indices values of 11% & 16% and 30% & 19%, respectively. In the case of NSE monthly TP load parameter interactions are relatively dominant compared to t otal TP load. 12 A , B and C . It was found that D_Sl_ON_a is the most important paramet er for all outputs (NSE monthly avgerage TN concentration, NSE monthly TN load, and t otal TN load) with corresponding ST i s ~ 73%, 71%, & 73%, and S i s ~ 46%, 21%, & 65%, respectively. The other important parameter, for a ll three outputs, was D_Sl_ON_b. Its S i s varied from <0.5% to 16% while corresponding ST i s varied from 24% to 45% between three outputs . The c ontributions of the remaining parameters to output variances were generally very small (< 5%) with notable exceptions of D_GW_vel and D_Sl_ON_Cb (dryland soluble organic nitrogen background concentration) in the case of NSE monthly TN load . Uncertainty Analysis PDFs, CDFs, 95% confide nce intervals (2.5% 97.5%) , and medians of 13, 4 14 and 4 - PAGE 125 125 15. None of the distributions were symmetric. Rather , most of them were considerably left skewed with the notable exception of t otal TN load ( F ig 4 15 C ) which was right skewed . 95% confide nce intervals for NSE daily flow, NSE monthly flow, and t otal flow volume were 0.26 0 .79, 0.869 0.907, and 3.96x10 8 4.24x10 8 (m 3 ), respectively. Corresponding medians were 0.55, 0 .898, and 4.17x10 8 (m 3 ) . Note that NSE daily flow had a tri modal PDF with modes at 0.31, 0.5, and 0.72 . Medians i.e. 50 percentile values for NSE monthly avg. TP con., NSE monthly TP load, and t otal TP load were ~ 0.5, 0.86, and 2.6x10 5 kg i.e. 260 metric tons (mt) . 95% confide nce intervals for these outputs were 2.1 0, 0.6 0.9, 1 25 3 50 mt . For NSE monthly avg. TN con., NSE monthly TN load, and t otal TN load median val ues were ~ 8, 0.8, and 830 mt, respectively. Confide nce intervals for these outputs were ~66% to 2%, 0.63 0.87, and 570 1340 mt. Note that t otal TN PDF load was found to be sk ewed to the right . Monte Carlo Filtering Earlier we noted that MCF was perfo r med considering multiple objectives and for this purpose addi tional outputs (likelihood function PBIAS outputs ) were also considered. It can be observed from the UA results that a ll NSE monthly avg. TP con. and NSE monthly avg. TN con. were less than the corresponding criteria (Table 4 3). For multi objective MCF these outputs and hence corresponding objectives were not included. This MCF exercise resulted in indetification of 318 parameter sets for which WAM/BLASROUTE passed all criteria. Posterior output distributions and corresponding ranges , except NSE monthly avg. TP and TN con., are presented in Figure s 4 16 A to G and Table 4 6 , respectively . For NSE daily flow, NSE monthly flow, t otal flow volume, NSE monthly TP load, and NSE monthly TN load the upper limit of the posterior distributions were the same as the prior distributions. The l ower limits of NSE daily flow and PAGE 126 126 NSE monthly flow correspon ded to the fil t ering criteria. Surprisingly, t posterior distribution range is similar to the prior MCF range. In the case of t otal TP load and t otal TN load , prior and posterior distributions are considerably different indicating reducti on in output uncertainty of these two outputs along with NSE outputs. Posterior parameter distributions are shown in Figures 4 17 to 4 19. Differences between prior and posterior distributions were tested by Kolgomorov Smirnov test, as explained in Chapte r 1, at 5% significance level. However, results indicated (not presented) that prior and posterior distributions are statistically different. Table 4 7 summarizes pr ior and post erior ranges of They do not differ much with notable exceptions of St_Sl_P_a and D_RO_kDF. For these two parameters posterior ranges are almost 50% of their respective prior ranges. Discussion Implications for WAM Applications, Pote n tial Model and Documentation Improvements Nutrient a ttenuation Nutrient a ttenuation m odel . Generalized attenuation equations for water quality constituents (Eqs 4 1 to 4 3) indicate that nutrient attenuation processes in the BLASROUTE module are independent of each other. In other words, BLASROUTE does no t account for biogeochemical processes once water quality constituents leave source cells. This is evident from EE ana lysis results as well. For example , for TP related outputs only P related parameters were identified as important. In fact EEs for the non P related water quality attenuation parameters (N, BOD, and TSS) were zero for P outputs. The s ame was true for TN, BOD, and TSS. F urther showed that all TP and TN outputs are almost completely dependent on a very few nutrient attenuation parameters, three (St_Sl_P_a, D_Sl_P_a, and D_Sl_P_b) in the case of TP and two (D_Sl_ON_a and D_Sl_ON_b) in the case of TN. It can be observed that d egree of the additivity PAGE 127 127 of TN and TP outputs varied substantially ( S i ~ 40% to 90%). Wher ever any TP and TN output was considerably non additive (say S i <70%) it was due to involvement of the above mentioned parameters in interactions. An argument could be that identification of these parameters as a the wide ranges assigned to them. However, large variations of their values was necessary due to their empirical nature (SWET, 2014). From the MCF exerci s e results it was noted that posterior parameter ranges were not much different from prior ranges provi ding evidence for the presence of equifinality (Beven and Binley, 1992) with respect to these parameters. T ypical values a b for the Floridian conditions are well established due to the long history of WAM applications in the state. However, for WAM a b attenuation parameters need t o be calibrated very carefully. Domina nt l and u se and r each t ype . WAM treats all land uses other than wetlands as dryland for the water quality attenuation, while for flow routing all land uses other than wetlands and depressions are treated as drylands. There were no depressions in the S 191 watershed. This implies that ~ 90% of the land u se was of the type dryland for BLASROUTE. Also, 83 out of 85 the In the EE parameter screening experiment 16 dryland related nutrient attenuation parameters were iden tified as important ( 12 very important + 4 less important) compared to 11 wetland related parameters (3 very important + 8 less important). The wetland related important parameters belonged to only two types freshwater marshes (3/6, very important/less i mpo r tant) and wetland prairies (0/2, very important/less impotant). These two wetlands types are the first and the second most abundant wetlands in S 191 with ~ 4.4% and 2% of land coverages, PAGE 128 128 respectively . showed that the nutrient attenuation parameters which can explain most of the output variance belong only to dryland. Among the time delay parameters, three parameters related to dryland were found to be important (3 very important and 0 less important) while o nly one parameter related to wetlands (wetland runoff constant delay W_RO_kDF ) was found to be important (one of the less important) in EE analysis . only 2 parameters (D_RO_kDF and D_GW_vel) had substantial influence o n model outputs, both o f them being related to dryland . Similar to dryland type parameters being identified as important, parameters belonging to only stream reaches were important. Among nitrogen species parameters related to organic nitrogen controlled TN outputs. From measured data it was found that concentrations of organic nitrogen were the highest among the various nitrogen species for S 191. Based on these observations it can be stated that the number of important parameters for a given land use or stream type is proportional to its coverage. In other words, for future WAM applications it would be most logical to tune parameters belonging to the dominant land use, reach type, and dominant nutrient specie (e.g. organic nitrogen in the case of TN) for model calibration. Model performance From uncertainty analysis results (Figures 4 13 A and B ) it was clear that WAM captures the dynamics of daily and monthly flow quite well as NSEs corresponding to them were greater than 0.65 for 40% and 100% of the sim ulation s , respectively. o z Carpena (2013) NSE = 0.65 is the threshold between acceptable and unacceptable results while according to Moriasi et al. (2007) NSE > 0.65 indicates good and very good model simulations. Note that for MCF , NS E > 0.75 were termed as behavio ral which is the threshold to identify very good model runs as per Moriasi et al . (2007). PBIAS values for both daily and month ly flow were in PAGE 129 129 the range 4.4% to + 5.8% indicating very good model performance s based on Moria si et al. (2007) according to which abs(PBIAS)<10% corresponds to very good model performance for flow outputs . NSEs for the monthly avg. TP and TN concentrations were negative f or over 98% of the simulations (Figures 4 14 and 4 15) . However, NSEs for mont hly TP and TN load outputs were much better with over 90% of the simulations having NSE > 0.65. This trend of improved model performance by switching from concentration to load at the monthly scale was also evident in daily outputs (results not presented). The r an ce in simulating nutrient concentrations but ability to capture load dynamics pretty well was not investigated further. However, the potential reasons for u nacceptable concentration but good load results are the unav ailability of continuous (daily) water quality data and the method o logy used to estimate measured and calculated loads . Note that the number of measured data points over the simulation period (Jan 2002 to Dec 2004) for P were 203 while for N were 51 . Correspondingly TP concentration results are relatively better than TN concentration results. The l oad calculation methodology used in this study is as per the procedure recommend ed by SFWMD (SWET personal communication). Note that N SE (eq. 4 8 ) is one of the most widely used goodness of fit measure / model performance indicator in H/WQ modeling . During recent decade s a number of authors have studied various aspects of NSE as a reliable model performance indicator, and have suggested modifications, and al ternatives for NSE (e.g. Legates and McCabe, 1999; McCuen et al., 2006; Jain and Sudheer, 2008; Gupta and King, 2011 etc.). One of the suggestions to improve the meaning of NSE as performance measure is to use only unique model observation data pairs (e.g. oz Carpena , 201 3 ). The s ame idea can be extended to other goodness of fit PAGE 130 130 measures. However, in this study we calculated NSE and PBIAS in a traditional manner i.e. considering all model observation data pairs. Secondly, measurement uncertai nty was not accounted while calculating goodness of fit measures. The framework for the incorporation of measurement uncertainty in model evaluation exists (e.g. Harmel and Smith, 2007; Harmel et al. oz Carpena, 2013). However, no study which simultaneously considered parametric and measurement uncertainty for model evaluation was found. NSE daily flow and D_RO_kDF The NSE daily flow distribution showed three distinct modes (Figure 4 13 A ). It may be speculated that the correspondin g output response surface has large plateaus which may pose a problem for some parameter optimization algorithms . This warranted an investigation. NSE daily fou nd that the tri modal behavior of NSE daily flow can be explained with only one parameter D_RO_kDF. A s catterplot of NSE daily flow against D_RO_kDF (Figure 4 20) indicates different behavior of WAM/BLASROUTE in four distinct parameter ranges 0 0.5, 0.5 1.5, 1.5 2.5, 2.5 3. The first mode (0.32) corresponded to D_RO_kDF = 2.5 3, the second mode at 0.5 to 1.5 2.5, while the last mode (0.72) corresponds to 0 1.5. This also explains why the posterior D_RO_kDF distribution range was much different from the prior MCF range. Looking at the D_RO_kDF ranges correspondin g to NSE daily values it may be speculated that D_RO_kDF values are rounded off by the WAM code. The s ame may be true for W_RO_kDF (wetland to runoff delay factor). However, due to the very small percentage of wetland coverage in this watershed , the effect may not be visible. I t is necessary to investigate t h is distin c t nature of WAM performance in the case of daily flow dynamics bu t such study is beyond the scope of this work. PAGE 131 131 Inconsistencies in WAM d ocumentation SWET Inc. (2014) has updated documentation of WAM to bring clarity in modeling processes as recommended by the peer review committee (Graham et al., 2009). However, during this model evaluation research we came across a few inconsiste ncies , some of which directly related to BLASROUTE parameters and input files are as mentione d below. 1. It was interesting to note that wetland to groundwater related time delay parameters did not affect routing and attenuation results at all, i.e. for these two parameters all EEs were zero. Careful inspection of flow distance grids indicated that there is no such grid for wetland to groundwater flow indicating that there is no groundwater flow from wetlands. The WAM manual does not explain thi s inconsistency i.e. presence of wetland to groundwater delay parameters and absence of corresponding distance grid. In other words W_GW_kDF and W_GW_vel are not at all parameters, but WAM manual does not explain this well. 2. For accurate reading of time del ay parameters by the WAM simulation code, delay parameters corresponding to dryland need to be specified on the last line of file. The s ame is true for the attenuation parameter file. WAM manual does not clear l y mention this. 3. Similar to the issue related to wetland to groundwater flow there is also an issue of time delay parameters related to depressions. For depressions there is only one distance grid dep ression to groundwater. Hence, in WAM there is provision of two delay parameters for depressions. These parameters are not explained explicitly in the WAM manual. For this study these parameters were not relevant as there were no depressions. 4. The n utrient attenuation equation presented in WAM manual (SWET, 2014) is in correct . The c orrect form of the equation (eq. 4 1) was determined after literature survey and discussion with the WAM developers. Also in the WAM documentation only one generalized att enuation equation is presented . This equation is cor rect for overland nutrient attenuation but incorrectly represents groundwater attenuation. The c orrect form of groundwater attenuation is presented in eq. 4 2 and was determined after discussion with WAM developers. Evaluation Techniques Comparison of EE a nd results SA of high dimensional models becomes a computationally expensive affair. The two step procedure for global SA used in this study was adopted to alleviate the burden on computational and time resources. Though this approach is not n ovel t o H/WQ studies , the PAGE 132 132 inclusion of more than 100 uncertain parameters in global SA is still rare. S ome recent studies have evaluated the efficiency of EE analysis to explore its potential to analysis (e.g. Herman et al., 2013 ; Wainwright et al., 2014 ). Similar to Chapter 3 , these studies used the same number of model parameters for both EE screening Though analysi s between ranks of the important parameters (Table 4 4 and 4 5) from EE analysis with 8). Overa ll correlations are high for all outputs. ble 4 4 and 4 5) are output. The o nly exceptions are NSE monthly TP load and NSE monthly TN load. However, were the top four On the other hand in the case of NSE monthly TN load , the model response is highly non additive ( S i ~ 30%) . Misclassificat ion of p arameters The original development of the EE method (Morris, 1991) was for qualitative separation of model parameters into important and non important classes. Even though the method has gone under subsequent modification for sensitivity measures (Campolongo et al., 2007, Campolongo et al., 2011) as well as sampling (Campolongo et al., 2007; Ruano et al., 2012; Chapter 3) the criteria for parameter segregation has not changed. In other words there are no established thresholds to classify parameter s into very important, somewhat important, and completely irrelevant categories. Identification of parameters as important or irrelevant is done by visual inspection of the separation of a parameter from the vertical axis of the * plot. Hence, modelers should be careful while interpreting EE based results while selecting parameters for optimization or quantitative SA i.e. classifying parameters into different categories . PAGE 133 133 In this study we observed some misclassification of paramete rs. For example for NSE daily flow n and dep_err were well separated from the origin in * indices indicated that they are very much irrelevant. The s ame was the case with n corresponding to NSE monthly TP load, St_Sl_ON_a corres ponding to NSE monthly avg. TN con., D_Sd_ON_a and FM_Sl_ON_Cb corres ponding to NSE monthly TN load . The only instance analysis was with RO_HYD. However, this inconsis tency RO_HYD is likely due to the different number of cate gories used for the two analys e s (4 vs. 20 hydrograph profiles, Figues C 1 and C 3) prompting the need to use a cautio u s approach while dealing with categorical par ameters. Identification of p arameter i nteraction and n on monotonic e ffects It is a standard practice to look at * alone to draw conclusions about the importance of parameter(s) on model output(s). * is always positive and hence cannot provide any information a change in a given parameter. Such information is useful to verify if the model code is correct or not. can provide this infor mation, though it suffers at times from inaccurate estimation of a importance . Morris (1991) suggested to plot two lines corresponding to = +/ 2SEM = +/ 2 /sqrt( r ). All the parameters that are outside the wedge formed by these two lin es on space can be considered to have significantly different from 0 and the corresponding parameter can be considered to be important by itself. Wainwright et al. (2014) took this further and tried to interpret a spac e in terms of non linearities and interactions. In the case of * space there is no published criteria, such as +/ 2SEM, that can give insights about relative contr ibutions of interaction and/ or non monotonicities and direct effects. In other words PAGE 134 134 it is important to look at both * and spaces instead of one to better understand EE results. Accordingly, we have plotted our EE results in both spaces (Figures 4 5 to 4 9). response to a parameter is monotonic over a defined distribution then for that parameter * = abs( ) should hold true. If that is not the case then it implies involvement of that particular parameter in a non monotonic response and/or parameter interactio ns. In figures 4 5 to 4 9 we have indicated the parameters that produce a non monotonic model response and/or are involved in interaction effects with the blue color while parameters with only a monotonic effect are indicated in black. It can be seen that all the parameters in black were outside the wedge formed by +/ 2SEM lines which is consistent with what Morris (1991) had proposed. Table 4 9 compares the ratio o f interaction effects ( ST i S i ) to ST i with EE based classification of parameters into blue a nd black. Note that we have considered parameters with ST i > 20%. It can be observed that for parameters involved in non monotonic/interaction effects , i.e. blue type , SI / ST ratio is generally high (>0.45) and for black type parameters its generally low (< 0.45) with the notable exception of D_Sl_ON_b for NSE monthly avg. TN con. Based on these results the black blue classification shows a good promise but need s further investigation . Normalization of s ensitivity m easures One of the major drawbacks of curre ntly used sensitivity measures in EE analysis is that the results for different outputs cannot be compared directly. Also, importance c an only be assessed relatively. One of the potential ways to overcome both these issues wo uld be to normalize EE based sensitivity measures. Summary In this work we have successfully performed sensitivity analysis of the routing and attenuation module of WAM treating all corresponding parameters as uncertain. This was the first ever sensitivity analysis study of WAM a model which is used very extensively in the PAGE 135 135 state of Florida in a number of important projects e.g. the Lake Okeechobee. The f indings of this study indicated the dependence of WAM on several emp irically derived water quality attenuation parameters. It was also found that model outputs were strongly dependent on the parameters of the most abundant land use in the study area. Uncertainty analysis and MCF manifested prospective issues with parameter optimization. Apart from model evaluation we identified potential avenues for improving analysis techniques especially the method of EE . It was found that presenting EE results in both * and spaces c an give insights about the changes in parameters (positive or negative) and involvement of a parameter in non mono tonic and/or interactive respons es . One major limitation of this study is that evaluation of WAM was performed only on BLASROUTE. However, such a simplified approach was necessary considering (a) the extremely large dimensionality of WAM, (b) absence of any formal SA/UA study of WAM, and (c) pronounced empiricism in the BLASROUTE module. Though useful, findings of this study do not po rtray holisti c picture about WAM and hence should be used with cauti on in future applications . PAGE 136 136 Table 4 1. Percentage land use distribution for the S 191 watershed Sr . Type % 1 Urban and built up 5.56 2 Pastures (all types) 60.74 3 Scrub and brushland 6.89 4 Wetlands 8.87 5 Agriculture 9.00 6 Hardwood and Conifer Mix 3.83 7 Others 5.11 PAGE 137 137 Table 4 2. Uncertainty analysis, and Monte Carlo Filtering Output EE Analysis Sobol' Analysis Uncertainty Analysis Monte Carlo Filtering NSE 1 daily flow YES YES YES YES NSE monthly flow YES YES YES YES Total flow volume YES YES YES NSE monthly avg. TP con. YES YES YES NSE monthly TP load YES YES YES YES Total TP load YES YES YES NSE monthly avg. TN con. YES YES YES NSE monthly TN load YES YES YES YES Total TN load YES YES YES Total BOD load YES Total TSS load YES PBIAS 2 daily flow YES PBIAS monthly flow YES PBIAS monthly avg. TP con. YES PBIAS monthly TP load YES PBIAS monthly avg. TN con. YES PBIAS monthly TN load YES 1 NSE Nash Sutcliff e Efficiency 2 PBIAS Pe r centage Bias PAGE 138 138 Table 4 3. Criteria 1 used for Monte Carlo Filtering (MCF) Flow Outputs TP Outputs TN Outputs NSE 2 >0.75 >0.75 > 0.75 PBIAS 3 = < +/ 10% = < +/ 25% = < +/ 25% 1 Based on Moriasi et al. (2007) 2 NSE Nash Sutcliff e Efficiency 3 PBIAS Percentage Bias PAGE 139 139 Table 4 4 . Summary of important WAM/BLASROUTE parameters for flow, BOD, and TSS related outputs . Numbers in parenthesis indicate parameter rank based on * . Please refer Table C 1 for parameter definitions. NSE_d aily_ flow NSE_m onthly_ flow Total flow volume Tot al TSS load 1 T ot al BOD load 1 Very Important Parameters D_RO_kDF (1) D_RO_kDF (1) D_GW_vel (1) D_TSS_Cb (1) D_Sl_BOD_Cb (1) n (2) n (2) dep_err (2) D_TSS_a (2) D_Sl_BOD_a (2) dep_err (3) D_GW_vel (3) FM_TSS_Cb (3) D_Sl_BOD_b (3) dep_err (4) D_TSS_b (4) FM_Sl_BOD_Cb (4) Less Important Parameters W_RO_kDF (5) wid_err (3) FM_TSS_a (5) FM_Sl_BOD_a (5) D_RO_vel (6) n (4) WP_TSS_Cb (6) WP_S l _ BOD_Cb (6) RO_HYD (7) FM_TSS_b (7) FM_Sl_BOD_b (7) corresponding parameters were included PAGE 140 140 Table 4 5 . Summary of important WAM/BLASROUTE parameters for TP and TN related outputs . Numbers in pare nthesis indicate parameter rank based on * . Please refer Table C 1 for parameter definitions. NSE_ monthly avg. TP con. NSE_ monthly TP load Total TP load NSE_ monthly avg. TN con. NSE_ monthly TN load Total TN load Very Important Parameters St_Sl_P_a (1) D_Sl_P_a (1) D_Sl_P_a (1) D_Sl_ON_a (1) D_Sl_ON_a (1) D_Sl_ON_a (1) D_Sl_P_a (2) St_Sl_P_a (2) St_Sl_P_a (2) D_Sl_ON_b (2) D_Sl_ON_b (2) D_Sl_ON_b (2) D_Sl_P_b (3) D_Sl_P_b (3) D_Sl_P_b (3) St_Sl_ON_a (3) D_GW_vel (3) St_Sl_P_Cb (4) n (4) D_Sd_ON_a (4) FM_Sl_ON_Cb (5) D_Sl_ON_Cb (6) Less Important Parameters D_GW_vel (5) dep_err (5) n (4) St_Sl_NO3_a (4) D_Sl_NO3_a (7) D_Sl_NO3_a (3) D_Sl_P_Cb (6) D_RO_kDF (6) D_Sl_P_Cb (5) dep_err (5) D_RO_kDF (8) FM_Sl_ON_Cb (4) dep_err (7) D_Sl_P_Cb (7) St_Sl_P_Cb (6) D_GW _ ON_a (6) n (9) D_Sl_ON_Cb (5) St_Sd_P_a (8) D_Sd_P_Cb (8) D_Sd_P_Cb (7) D_GW_vel (7) FM_Sl_ON_a (10) D_GW_vel (6) D_Sd_P_Cb (9) St_Sl_P_Cb (9) D_Sl_NO3_a (8) St_Sl_ON_a (11) St_Sd_P_Cb (10) D_Sd_P_a (10) D_Sd_P_a (11) PAGE 141 141 Table 4 6. Summary of WAM/BLASROUTE outputs posterior to Monte Carlo Filtering (MCF) Output Minimum Maximum NSE daily flow 0.750 0.802 NSE monthly flow 0.886 0.908 Total flow volume (m 3 ) 3.9E+08 4.3E+08 NSE monthly TP load 0.832 0.894 Total TP load (kg) 2.3E+05 3.6E+05 NSE monthly TN load 0.750 0.890 Total TN load (kg) 6.2E+05 8.6E+05 PAGE 142 142 Table 4 7. Summary of parameter ranges pre and post Monte Carlo Filtering (MCF) Parameter Prior Po sterior Minimum Maximum Minimum Maximum St_Sl_NO3_a 0.0005 0.05 0.000585 0.049575 St_Sl_ON_a 0.0001 0.01 0.000103 0.009710 St_Sl_P_a 0.004 0.4 0.004198 0.200931 St_Sl_P_Cb 0.02 0.2 0.020684 0.196757 St_Sd_P_a 0.002 0.2 0.002104 0.188238 St_Sd_P_Cb 0.001 0.1 0.001005 0.078196 FM_Sl_ON_a 0.00001 0.001 0.000010 0.000795 FM_Sl_ON_Cb 0.5 3 0.560995 2.996546 D_Sl_NO3_a 0.00001 0.001 0.000011 0.000998 D_Sl_ON_a 0.00001 0.001 0.000032 0.000970 D_Sl_ON_b 0.45 0.85 0.455372 0.846674 D_Sl_ON_Cb 0.05 0.5 0.052568 0.493828 D_Sl_P_a 0.000001 0.0001 0.000001 0.000096 D_Sl_P_b 0.45 0.85 0.450415 0.833582 D_Sl_P_Cb 0.005 0.3 0.005041 0.284662 D_Sd_ON_a 0.000001 0.0001 0.000001 0.000095 D_Sd_P_a 0.000001 0.0001 0.000001 0.000090 D_Sd_P_Cb 0.005 0.15 0.005612 0.149951 D_GW_ON_a 0.4 0.95 0.407321 0.949557 RO_HYD 1 20 1 20 D_RO_vel 1000 40000 1144 39224 D_RO_kDF 0 3 0.017120 1.430838 D_GW_vel 2 200 2.210720 196.805626 W_RO_kDF 0 3 0.006188 2.990912 dep_err 0.2 0.2 0.178584 0.197239 wid_err 0.2 0.2 0.199860 0.188676 n 0.02 0.08 0.020709 0.079738 PAGE 143 143 Table 4 8. Correlation coefficients between ranks of parameters obtained from EE screening and 1 . Output Very Important Parameters All Parameters NSE daily flow 0.982 0.982 NSE monthly flow 0.990 0.640 Total flow volume 1.000 1.000 NSE monthly avg. TP con. 1.000 0.776 NSE monthly TP load 0.800 0.872 Total TP load 1.000 0.817 NSE monthly avg. TN con. 0.982 0.643 NSE monthly TN load 0. 829 0.932 Total TN load 1.000 0.486 1 Very Important parameters include only those parameters which were identified as very important in EE analysis PAGE 144 144 Table 4 9 . Comparis 1 and proposed criteria for identification of non monotonic/interaction effect from EE anal ysis Output Parameter SI i / ST i 2 EE classification indicator 3 NSE daily flow D_RO_kDF 0.0558 Black NSE mon thly flow RO_HYD 0.7334 Blue D_RO_kDF 0.5179 Blue n 0.6247 Blue Total flow volume D_GW_vel 0.0276 Black NSE monthly avg. TP con . St_Sl_P_a 0.4883 Blue D_Sl_P_a 0.7622 Blue D_Sl_P_b 0.6259 Blue NSE monthly TP load St_Sl_P_a 0.3928 Blue D_Sl_P_a 0.4021 Black D_Sl_P_b 0.6224 Blue Tot al TP load St_Sl_P_a 0.2640 Black D_Sl_P_a 0.1526 Black NSE monthly avg. TN c on. D_Sl_ON_a 0.3487 Black D_Sl_ON_b 0.6062 Black NSE monthly TN load D_Sl_ON_a 0.7197 Blue D_Sl_ON_b 0.9989 Blue Tot al TN load D_Sl_ON_a 0.1154 Black D_Sl_ON_b 0.3339 Black 1 Only those parameters for which ST i 2 SI i = ST i S i 3 If for a parameter * parameter can be classified as not involved in non monotonic responses/interactions and is coded black else blue . PAGE 145 145 Figure 4 1. Determination of unique model cells in WAM (adopted fro m WAM Documentation, S WET 2014 ). PAGE 146 146 Figure 4 2. Schematic representation the BUCSHELL module ( adopted from WAM Documentation, SWET 2014 ). PAGE 147 147 Figure 4 3. Schematic of the BLASROUTE sub module flow directions and distance grids (adopted from Bottcher et al. 2012 ). PAGE 148 148 Figure 4 4. Maps showing A ) Location and topography of S 191 watershed, B ) Land use within S 191 watershed and C ) Distribution of soils within S 191 watershed . A B C PAGE 149 149 Figure 4 5. Elementary Effects sensitivity analysis results for the flow related outputs A ) NSE B C ) NSE monthly flow in D E ) Total flow volume in space, and F space . F E C D A B PAGE 150 150 Figure 4 6. Elementary Effects sensitivity analysis results for the TP related outputs A ) NSE monthly average TP concentration B ) NSE mon thly average TP concentration C ) NSE monthly TP load D ) NSE monthly TP load E ) Total TP load in F ) Total TP load in space. A B C D E F PAGE 151 151 Figure 4 7. Elementary Effects sensitivity analysis results for the TN rel ated outputs A ) NSE monthly average TN concentrations B ) NSE monthly average TN concentrations C D ) NSE monthly TN load in E ) Total TN load in F ) Total TN load in space . A B C D F E PAGE 152 152 Figure 4 8. Elementary Effects sensitivity analysis results for the BOD related outputs A ) Total B ) Total BOD load in A B PAGE 153 153 Figure 4 9. Elementary Effects sensit ivity analysis results for the TSS related outputs A ) Total B ) Total TSS load in A B PAGE 154 154 Figure 4 10. effect sensitivity indices ( S and ST ) , for flow related outputs A ) NSE daily flow, B ) NSE monthly flow, and C ) Total flow volume (m 3 ). Parameters for which both S and ST were <0.5% or 0.005 were not plotted. A B C PAGE 155 155 Figure 4 11. Results of ( S and ST ) , for TP related outputs A ) NSE monthly avg. TP con., B ) NSE monthly TP load, and C ) Total TP load (kg). Parameters for which both S and ST were <0.5% or 0.005 were not plotted. C B A PAGE 156 156 Figure 4 12. ( S and ST ) , for TN related outputs A ) NSE monthly avg. TN con., B ) NSE monthly TN load, and C ) Total TN load (kg). Parameters for which both S and ST wer e <0.5% or 0.005 were not plotted. C A B PAGE 157 157 Figure 4 1 3 . Probability Density Function (PDF), Cumulative Density Function (CDF), 95% confid e nce interval (2.5% 97.5%), and median values for WAM/BLASROUTE outputs A ) NSE daily flow, B ) NSE monthly flow, and C ) Total flow volume (m 3 ). C B A PAGE 158 158 Figure 4 1 4 . Probability Density Function (PDF), Cumulative Den sity Function (CDF), 95% confide nce interval (2.5% 97.5%), and median values for WAM/BLASROUTE outputs A ) NSE monthly avg. TP con., B ) NSE monthly TP load, and C ) Total TP load (kg). C A B PAGE 159 159 Figure 4 15. Probability Density Function (PDF), Cumulative Den sity Function (CDF), 95% confide nce interval (2.5% 97.5%), and median values for WAM/BLASROUTE outputs A ) NSE monthly avg. TN con., B ) NSE monthly TN load, and C ) Total TN load (kg). C B A PAGE 160 160 Figure 4 16. Posterior distributions of WAM/BLASROUTE outputs used in Monte Carlo Filtering A ) NSE daily flow, B ) NSE monthly flow, C ) Total flow volume (m 3 ), D ) NSE monthly TP load, E ) Total TP load (kg), F ) NSE monthly TN load, and G ) Total TN load (kg). C A E G F D B PAGE 161 161 Figure 4 17. sensitivity analysis and uncertainty analysis A ) St_Sl_NO3_a, B ) St_Sl_ON_a, C ) St_Sl_P_a, D ) St_Sl_P_Cb, E ) St_Sd_P_a, F ) St_Sd_P_Cb, G ) FM_Sl_ON_a, H ) FM_Sl_ON_Cb, I ) D_Sl_NO3_a, J ) D_Sl_ON_a, K ) D_Sl_ON_b, L ) D_Sl_ON_Cb A B C D E F G H I J K L PAGE 162 162 Figure 4 18. sensitivity analysis and uncertainty analysis A ) D_Sl_P _a, B ) D_Sl_P_b , C ) D_Sl_P_Cb , D ) D_Sd_ON_a , E ) D_Sd _P_a, F ) D _Sd_P_Cb, G ) D_GW_ON_a , H ) RO_HYD , I ) D_RO_vel, J ) D_ RO_kDF , K ) D_ GW_vel , L ) W_RO_kDF A B C D E F G H I J K L PAGE 163 163 Figure 4 19. sensitivity analysis and uncertainty analysis A ) dep_err, B ) wid_err, C ) n A B C PAGE 164 164 Figure 4 20. Scatterplot of NSE daily flow and D_RO_kDF PAGE 165 165 CHAPTER 5 CONCLUSIONS The use of hydrologic and water quality models (H/WQ) as regulatory decision making tools has increased substantially over the last few decades. For accurate and well informed decision making it is necessary to test H/WQ models for their reliability. Regulatory agencies such as United State s Environmental Protection Agency (EPA, 2009) and European Commission (EC, 2009) have emphasized the inclusion of risk assessment in environmental model applications . Uncertainty and sensitivity analyses (UA/SA) are important model evaluation tools that are used for risk reduction . In spite of their importance in model building UA/SA have not been historically imp lemented by the H/WQ community. Though th is has changed during the last ~ 15 years , there a re a number of aspects of UA/SA that are not yet well understood. This research was centered on the theme of understanding, improving, and applying SA techniques to H /WQ models. The literature review in C ha pter 1 indicated a number of problems r elated to SA/UA techniques that are generally applied in H/WQ modeling applications. Two specific issues, typical to SA in H/WQ model evaluations , include (a) parameter correlations and (b) high dimensionality. Most of the currently available SA/UA techniques make the assumption that model parameters are orthogonal (not correlated). Though there are some less elegant SA techniques to handle such situation s , para meter correlation is often ignored. The h igh dimensionality of large scale H/WQ models incre a ses the computational demand of SA application s beyond the point of practical implem en tation. To address the is s ue of high dimensionality t here are ways to qualita tively identify important parameters through screening processes. However, the most widely accepted and used technique for this purpose the method PAGE 166 166 of Elementary Effects/ Morris method needs improvement. We have addresse d these issues in this study by div iding them into three objectives: To perform global sensitivity and uncertainty analysis of a simple drought index model with correlated parameters using two different sensitivity analysis techniques . To develop and evaluate a new multi criteria trajectory based parameter sampling strategy for the screening type sensitivity analysis method of Elementary Effects . To implement the newly developed parameter sampling strategy as a part of a sensitivity and uncertainty evaluation framework for a high dimensional watershed scale model . Global sensitivity and uncertainty analysis of a simple drought index model A gricultural R eference I ndex for D rought (ARID) was performed successfully in Chapter 2 considering four textural soil classes and five locations across the southeast USA . One of the challenges was the presence of correlated model parameters when employing a variance based sensitivity technique. To account for this issue an innovative strategy was applied to combine the two corre lated parameters into another surrogate parameter. Results of partial correlation analysis (less elegent SA technique) and the (variance based technique) were found to be comparable . SA results (both techniques) indicate d that ARID is a mo notonic and linear model and is highly influenced by root zone depth followe d by soil hydraulic properties. These results have important implications while forecasting agricul tural water demands using ARID (e.g. AgroClimate tools). The use of appropriate e valuation methodologies is of utmost importance when H/WQ models are used to support decision making and the selection of a suitable parameter sampling strategy is a multi faceted problem. We develop e d a new multi criteria trajectory based parameter sampli ng strategy (SU) for calculating sensitivity measures based on the method of EEs for the parameter screening exercises (Chapter 3) . The main concept around which the SU is built is honoring the theoretical parameter distribution (i.e. discrete uniform distribution) along PAGE 167 167 with maintaining good spread of parameter samples. While SU outperform ed some of the trajectory based benchmark sa mpling strategies in most of the evaluation criteria its implementation is flexible to control trajectory spread by selecting the desired value of the oversampling size. The advantage of SU as a sampling strategy is founded on combined sampling criteria. T his work on improving the EE analysis technique is of paramount significance in the context of reliable screening of important model parameters for reducing the dimensionality of H/WQ model s in the calibration validation process. In Chapter 4 we implemente d SU (sampling s trategy developed in Chapter 3) as a part of a SA framework to analy ze the Watershed Assessment Model (WAM) for parametric sensitivity and uncertainty . It is a fully distributed H/WQ model with a large number of input parameters. WAM is one of the most widely used watershed model s for environmental policy making and evaluation (e.g. Total Maximum Daily Loads formualtion) in southern Florida a region of complex hydro ecological systems. This was the first ever global SA/UA evaluation study of WAM. The test case watershed, S 191 basin in south Florida, is one of the major contributors of phosphorus to Lake Okeechobee and has been previously modeled using WAM. The number of WAM/BLASROUTE parameters treated uncertain in this study were ~ 300. I t was found that the hydrologic and water quality routing algorithms in WAM are highly dependent on conceptual equations/ empirical parameters . The depiction of biogeochemistry of water quality constituents is highly simplified and routing/attenuation of an individual constituent is independent of others. Further, the presence of equifinality was evident. Findings of this study will help WAM user s and developers to take necessary measures to improve its applications in environmental policy formulation and evaluation. PAGE 168 168 Limitations and future r esearch . The three specific objectives of this study, mentioned previously, were successfully achi e ved. Howe ver, the corresponding individual studies were not free of assumptions and limitations. They are discussed in detail in the individual chapters. The m ain l i mitations and scopes of potential future research for the individual chapters are as follows Chapter 2 Climat e related inputs to ARID were not treated as uncertain for a given location. The e ffect of meteorology on ARID was studied implicitly by performing analysis for five different locations. Though the approach to handle parameter correlation for this study is applicable for all types of models, more elegant variance based SA/UA techniqu es should be developed that can explicitly account for parameter non orthogonality. Chapter 3 Only t rajectory based sampling strategies for the EE method were evaluated and compared . SU combines uniformity and spread of the generated samples as the sampling criteria. In the future other sampling combinations e.g. Man hattan distance criteria instead of E uclidean distance should be explored. Chapter 4 Only the flow and nutrient routing module of WAM i.e. BLASROUTE was evaluated for parametric sensitivity and uncertainty. Considering this , findings of this study should be used carefully . For the holistic pi cture about WAM uncertainty , future efforts (if possible) should simultaneously treat all WAM parameters and input data (meteorology, land use, soil etc.) as uncertain. In spite of its usefulness as a reliable parameter screening technique at times the method of EE requires personal judgment by the modeler to classify parameters into important and irrelevant categories. Future efforts to improve the method of EE, which can help reducing the computational burden on variance based global SA/UA, should be d irected towards standardization of screening criteria along with normalization of sensitivity measures. PAGE 169 169 This research focused on model evaluation through uncertainty and sensitivity analysis of hydrologic and water quality models will be continued and impr oved by addressing shortcomings of this work as mentioned in above. PAGE 170 170 APPENDIX A ALGORITHM AND MATLAB PROGRAM FOR THE SAMPLING FOR UNIFORMITY SAMPLING STRATEGY Matlab package to generate parameter trajectories for the method of elementary effects using the new strategy developed in this work is available for free download from http://abe.ufl.edu/carpena/software/SUMorris.shtml . Below is the main matlab function that needs to be called by the user from the command line in order to generate parameter sample. This code further calls a number of other matlab functions (not included here). It should be noted that in this program numbe r of levels are hard coded to 4. function [Traj,D,time] = SU_Sampling(NumFact,Q,r,trunc) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % SU_Sampling function generates parameter samples to perform parameter screening using the % method of Elementary Effe cts % This code is based on logic Sampling for Uniformity developed by Khare et al. % There are three inputs to this function: % NumFact = Number of Input Factor or Parameters to b e Sampled % Q = Oversampling Size % r = Number of Trajectories (should be an even number) % trunc = truncation to be done or not: 0 no truncation, 1 truncation at 0.125 and 0.875 % p = 4; Number of levels (Hard coded) % % Function outputs consist of: % (1) Traj: sample in th e form of 'r' traje ctories sampled from unit % hyperspace truncated at 0.125 and 0 .875 % (2) D: Euclidean Distance between traject ories, a measure of sample % spread % (3) time: time required for sample genera tion in seconds %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% c = unifrnd(1,10,1,1);%To induce randomness in trajectory generation c = ceil(c); if mod(r,2)==0 %Matrix init ilization %D(1,1)=zeros; D1(1:Q,1)=zeros; D2(1:Q,1)=zeros; %R(1:1,1)=zeros; time(1,1)=zeros; traj1(1:(NumFact+1)*r,1:NumFact)=zeros; PAGE 171 171 traj2(1:(NumFact+1)*r,1:NumFact)=zeros; Sampling_Stat(c,2); %Used to induce randomness in trajectory generation tic; %start of sample generation % for loop below repeats the sample generation experiment 'N' times and % selects the trajectories with highest spread for ii=1:Q %Function Sampling Stat generates 'r' trajectories for 'NumFact' %parameters and corresponding Euclidean distance 'D' [traj2,D2(ii,1)]=Sampling_Stat(NumFact,r); if ii==1 traj1=traj2; D1(ii,1)=D2(ii,1 ); else if D2(ii,1)>D1(ii,1) traj1=traj2; D1(ii,1)=D2(ii,1); end end end time=toc;% end of sample generation Traj=traj1; D=D1(Q,1); if trunc==0 [Traj_tr]=TR AJ_TRANS(Traj); Traj=Traj_tr; end R=Repeat_Count(Traj);%Check to find if there is any repetition of sample point if R(1,1)>0 fprintf(' \ n Caution: One or more sample points are repeated \ n'); end else % This program requir es 'r' to be an even number % If 'r' is not even, no sample is generated and following error is % flashed fprintf(' \ n Error: Number of Trajectoris (r) should be even \ n'); end PAGE 172 172 Figure A 1. Flowchart for Sampling for Uniformity (SU) to generate sample ( r trajectories) for k parameters with the resampling size of Q . PAGE 173 173 APPENDIX B TEST FUNCTIONS Details of test functions and the sets of important parameters are presented in this appendix. function and Annoni, 2010 etc.). It is a strongly non linear and non monoton ic function (JCR EC, 2010). Sets of important parameter, in the descending order of ranks, for the G function configurations in this study G15: T ST = { X 1 , X 2 , X 3 , X 4 }, #T = 4 G20: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 }, #T = 7 G35: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 }, #T = 10 G50: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 , X 11 }, #T = 11 G80: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 , X 11 , X 12 , X 13 }, #T = 13 G100: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 , X 11 , X 12 , X 13 , X 14 , X 15 , X 16 , X 17 , X 18 }, #T = 18 Modified G or G* function The G function was modified by Saltelli et al. (2010) to add complexity and flexibility (curvature, shift, etc.) to the function. The G* function has been used consistently in recently published works (e.g. C ampolongo et al., 2011). Sets of important parameter, in the descending order of ranks, for the G * function configurations in this study G*15: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 , X 11 , X 12 }, #T = 12 G*20: T ST = { X 2, X 8 , X 11 , X 18 }, #T = 4 G*35: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 }, #T = 10 PAGE 174 174 G*50: T ST = { X 1 , X 2 , X 4 , X 3 , X 6 , X 7 } , #T = 6 G*80: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 }, #T = 8 G*100: T ST = { X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 , X 9 , X 10 , X 11 }, #T = 11 B function The B function is a non additive test model developed by Saltelli et al. (2008). It has been used as a test case in a number of recent studies. Sets of important parameter, in the descending order of ranks, for the B function configurations in this study B20: T ST = { X 1 0 , X 20 , X 6 , X 9 , X 19 , X 5 , X 16 }, #T = 7 B50: T ST = { X 17 , X 19 , X 21 , X 23 , X 25 , X 50 , X 46 , X 16 , X 48 , X 44 , X 42 , X 18 , X 20 }, #T = 13 B80: T ST = { X 80 , X 40 , X 39 , X 79 , X 38 , X 78 , X 37 , X 77 , X 36 , X 76 , X 35 , X 75 , X 34 , X 74 }, #T = 14 B100: T ST = { X 100 , X 99 , X 98 , X 97 , X 96 , X 95 , X 94 }, #T = 7 M function The 20 parameter polynomial test function proposed by Morris (1991) has been the most regularly used test case in parameter screening experiments. The first and total effects sensitivity indices were calculated in an independent exper following parameters were identified as important T ST = { X 1 , X 4 , X 2 , X 9 , X 10 , X 8 , X 5 , X 3 , X 6 , X 7 }, #T = 10 O function The O function is a 15 parameter test model propose Model coefficient values can be downloaded from http://www.jeremy oakley.staff.shef.ac.uk/psa_example.txt . Sets of important parameters were identified from the analytical sensitivit y indices in Saltelli et al. (2008). T ST = { X 15 , X 11 , X 12 , X 13 , X 14 , X 9 , X 8 } , #T = 7 PAGE 175 175 Matlab (MathWorks Inc., 2012) programs for the first three functions (G, G* and B) were downloaded from JRC, EC ( http://ipsc.jrc.ec.europa.eu/?id=756 ) while the last two functions were coded by us for this study. Additional information about these test functions and codes can PAGE 176 176 APPENDIX C WAM/BLASROUTE PARAMETERS Table C 1. Summary of WAM/BLASROUTE parameters and their probabilistic distributions used for the purpose of parameter screening using the method of Elementary Effects Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 St_Sl_BOD_a Stream sol. BOD attn. 0.0001 0.00001 0.001 l og 10 u ni St_Sl_BOD_Cb ppm Stream sol. BOD bckgrnd. C on. 0.5 0.3 0.7 u ni St_Sl_NO3_a Stream sol. NO3 att n 0.005 0.0005 0.05 l og 10 u ni St_Sl_NO3_Cb ppm Stream sol. NO3 bckgrnd. C on. 0.02 0.001 0.02 l og 10 u ni St_Sl_NH4_a Stream sol. NH4 attn. 0.00001 0.000001 0.0001 l og 10 u ni St_Sl_NH4_Cb ppm Stream sol. NH4 bckgrnd. Con. 0.015 0.001 0.035 l og 10 u ni St_Sl_ON_a Stream sol. OrgN attn. 0.001 0.0001 0.01 l og 10 u ni St_Sl_ON_Cb ppm Stream sol. OrgN bckgrnd. Con. 0.2 0.15 1 u ni St_Sl_P_a Stream sol. P attn. 0.04 0.004 0.4 l og 10 u ni St_Sl_P_Cb ppm Stream sol. P bckgrnd. Con. 0.2 0.02 0.2 l og 10 u ni St_TSS_a Stream sus. solids attn. 0.00001 0.000001 0.0001 l og 10 u ni St_TSS_Cb ppm Stream sus. solids bckgrnd. Con. 4 2 6 u ni St_Sd_NH4_a Stream sed. NH4 attn. 0.00001 0.000001 0.0001 l og 10 u ni St_Sd_NH4_Cb ppm Stream sed. NH4 bckgrnd. Con. 0.007 0.001 0.01 u ni St_Sd_ON_a Stream sed. OrgN attn. 0.00001 0.000001 0.0001 l og 10 u ni St_Sd_ON_Cb ppm Stream sed. OrgN bckgrnd. Con. 0.07 0.05 0.1 u ni St_Sd_P_a Stream sed. P attn. 0.02 0.002 0.2 l og 10 u ni St_Sd_P_Cb ppm Stream sed. P bckgrnd. Con. 0.07 0.001 0.1 l og 10 u ni Cn_Sl_BOD_a Canal sol. BOD attn. 0.0001 0.00001 0.001 l og 10 u ni Cn_Sl_BOD_Cb ppm Canal sol. BOD bckgrnd. C on. 0.5 0.3 0.7 u ni Cn_Sl_NO3_a Canal sol. NO3 att n. 0.01 0.001 0.1 l og 10 u ni Cn_Sl_NO3_Cb ppm Canal sol. NO3 bckgrnd. con. 0.01 0.001 0.02 l og 10 u ni Cn_Sl_NH4_a Canal sol. NH4 attn. 0.00001 0.000001 0.0001 l og 10 u ni Cn_Sl_NH4_Cb ppm Canal sol. NH4 bckgrnd. Con. 0.015 0.001 0.035 l og 10 u ni Cn_Sl_ON_a Canal sol. OrgN attn. 0.001 0.0001 0.01 l og 10 u ni Cn_Sl_ON_Cb ppm Canal sol. OrgN bckgrnd. Con. 0.2 0.15 1 u ni Cn_Sl_P_a Canal sol. P attn. 0.0015 0.00015 0.015 l og 10 u ni Cn_Sl_P_Cb ppm Canal sol. P bckgrnd. Con. 0.02 0.02 0.2 l og 10 u ni Cn_TSS_a Canal sus. solids attn. 0.00001 0.000001 0.0001 l og 10 u ni Cn_TSS_Cb ppm Canal sus. solids bckgrnd. Con. 4 2 6 u ni Cn_Sd _ NH4_a Canal sed. NH4 attn. 0.00001 0.000001 0.0001 l og 10 u ni Cn_Sd _ NH4_Cb ppm Canal sed. NH4 bckgrnd. Con. 0.007 0.001 0.01 u ni Cn_Sd _ ON_a Canal sed. OrgN attn. 0.00001 0.000001 0.0001 l og 10 u ni Cn_Sd _ ON_Cb ppm Canal sed. OrgN bckgrnd. Con. 0.07 0.05 0.1 u ni Cn_Sd _ P_a Canal sed. P attn. 0.0015 0.00015 0.015 l og 10 u ni Cn_Sd _ P_Cb ppm Canal sed. P bckgrnd. Con. 0.05 0.001 0.1 l og 10 u ni PAGE 177 177 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 BS_Sl_BOD_a B.S. sol. BOD attn. 0.0003 0.00003 0.003 l og 10 u ni BS_Sl_BOD_b B.S. sol. BOD attn. 0.5 0.3 0.8 u ni BS_Sl_BOD_Cb ppm B.S. sol. BOD bckgrnd. Con. 0.5 1 5 u ni BS_Sl_NO3_a B.S. sol. NO3 attn. 0.0002 0.00001 0.001 l og 10 u ni BS_Sl_NO3_b B.S. sol. NO3 attn. 0.6 0.3 0.8 u ni BS_Sl_NO3_Cb ppm B.S. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni BS_Sl_NH4_a B.S. sol. NH4 attn. 0.0002 0.00001 0.001 l og 10 u ni BS_Sl_NH4_b B.S. sol. NH4 attn. 0.6 0.3 0.8 u ni BS_Sl_NH4_Cb ppm B.S. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni BS_Sl_ON_a B.S. sol. OrgN attn. 0.0002 0.00001 0.001 l og 10 u ni BS_Sl_ON_b B.S. sol. OrgN attn. 0.6 0.3 0.8 u ni BS_Sl_ON_Cb ppm B.S. sol. OrgN bckgrnd. Con. 1.6 0.5 3 u ni BS_Sl_P_a B.S. sol. P attn. 0.0000085 0.000001 0.0001 l og 10 u ni BS_Sl_P_b B.S. sol. P attn. 0.65 0.3 0.8 u ni BS_Sl_P_Cb ppm B.S. sol. P bckgrnd. Con. 0.1 0.005 0.3 l og 10 u ni BS_TSS_a B.S. sol. Sus. Solids attn. 0.0003 0.00003 0.003 l og 10 u ni BS_TSS_b B.S. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni BS_TSS_Cb ppm B.S. sol. Sus. Solids bckgrnd. Con. 10 2 15 u ni BS_Sd_NH4_a B.S. sed. NH4 attn. 0.0003 0.00003 0.003 l og 10 u ni BS_Sd_NH4_b B.S. sed. NH4 attn. 0.5 0.3 0.8 u ni BS_Sd_NH4_Cb ppm B.S. sed. NH4 bckgrnd. Con. 0.01 0.001 0.05 l og 10 u ni BS_Sd_ON_a B.S. sed. OrgN attn. 0.0003 0.00003 0.003 l og 10 u ni BS_Sd_ON_b B.S. sed. OrgN attn. 0.5 0.3 0.8 u ni BS_Sd_ON_Cb ppm B.S. sed. OrgN bckgrnd. Con. 0.1 0.015 0.15 l og 10 u ni BS_Sd_P_a B.S. sed. P attn. 0.00003 0.000003 0.0003 l og 10 u ni BS_Sd_P_b B.S. sed. P attn. 0.5 0.3 0.8 u ni BS_Sd_P_Cb ppm B.S. sed. P bckgrnd. Con. 0.1 0.005 0.15 l og 10 u ni BS_GW_NO3_a B.S. groundwater NO3 attn. 0.8 0.4 0.95 u ni BS_GW_NH4_a B.S. groundwater NH4 attn. 0.8 0.4 0.95 u ni BS_GW_ON_a B.S. groundwater OrgN attn. 0.8 0.4 0.95 u ni BS_GW_P_a B.S. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni SLS_Sl_BOD_a S.L.S. sol. BOD attn. 0.0003 0.00003 0.003 l og 10 u ni SLS_Sl_BOD_b S.L.S. sol. BOD attn. 0.5 0.3 0.8 u ni SLS_Sl_BOD_Cb ppm S.L.S. sol. BOD bckgrnd. Con. 0.5 1 5 u ni SLS_Sl_NO3_a S.L.S. sol. NO3 attn. 0.00023 0.00001 0.001 l og 10 u ni SLS_Sl_NO3_b S.L.S. sol. NO3 attn. 0.6 0.3 0.8 u ni SLS_Sl_NO3_Cb ppm S.L.S. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni SLS_Sl_NH4_a S.L.S. sol. NH4 attn. 0.00023 0.00001 0.001 l og 10 u ni SLS_Sl_NH4_b S.L.S. sol. NH4 attn. 0.6 0.3 0.8 u ni PAGE 178 178 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 SLS_Sl_NH4_Cb ppm S.L.S. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni SLS_Sl_ON_a S.L.S. sol. OrgN attn. 0.00023 0.00001 0.001 l og 10 u ni SLS_Sl_ON_b S.L.S. sol. OrgN attn. 0.6 0.3 0.8 u ni SLS_Sl_ON_Cb ppm S.L.S. sol. OrgN bckgrnd. Con. 1.6 0.5 3 u ni SLS_Sl_P_a S.L.S. sol. P attn. 0.000009 0.000001 0.0001 l og 10 u ni SLS_Sl_P_b S.L.S. sol. P attn. 0.65 0.3 0.8 u ni SLS_Sl_P_Cb ppm S.L.S. sol. P bckgrnd. Con. 0.05 0.005 0.3 l og 10 u ni SLS_TSS_a S.L.S. sol. Sus. Solids attn. 0.0003 0.00003 0.003 l og 10 u ni SLS_TSS_b S.L.S. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni SLS_TSS_Cb ppm S.L.S. sol. Sus. Solids bckgrnd. Con. 5 2 15 u ni SLS_Sd_NH4_a S.L.S. sed. NH4 attn. 0.0003 0.00003 0.003 l og 10 u ni SLS_Sd_NH4_b S.L.S. sed. NH4 attn. 0.5 0.3 0.8 u ni SLS_Sd_NH4_Cb ppm S.L.S. sed. NH4 bckgrnd. Con. 0.005 0.001 0.05 l og 10 u ni SLS_Sd_ON_a S.L.S. sed. OrgN attn. 0.0003 0.00003 0.003 l og 10 u ni SLS_Sd_ON_b S.L.S. sed. OrgN attn. 0.5 0.3 0.8 u ni SLS_Sd_ON_Cb ppm S.L.S. sed. OrgN bckgrnd. Con. 0.05 0.015 0.15 l og 10 u ni SLS_Sd_P_a S.L.S. sed. P attn. 0 0.000003 0.0003 l og 10 u ni SLS_Sd_P_b S.L.S. sed. P attn. 0.5 0.3 0.8 u ni SLS_Sd_P_Cb ppm S.L.S. sed. P bckgrnd. Con. 0.05 0.005 0.15 l og 10 u ni SLS_GW _ NO3_a S.L.S. groundwater NO3 attn. 0.8 0.4 0.95 u ni SLS_GW _ NH4_a S.L.S. groundwater NH4 attn. 0.8 0.4 0.95 u ni SLS_GW _ ON_a S.L.S. groundwater OrgN attn. 0.8 0.4 0.95 u ni SLS_GW _ P_a S.L.S. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni WH_Sl_BOD_a W.H. sol. BOD attn. 0.0003 0.00003 0.003 l og 10 u ni WH_Sl_BOD_b W.H. sol. BOD attn. 0.5 0.3 0.8 u ni WH_Sl_BOD_Cb ppm W.H. sol. BOD bckgrnd. Con. 0.5 1 5 u ni WH_Sl_NO3_a W.H. sol. NO3 attn. 0.00023 0.00001 0.001 l og 10 u ni WH_Sl_NO3_b W.H. sol. NO3 attn. 0.6 0.3 0.8 u ni WH_Sl_NO3_Cb ppm W.H. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni WH_Sl_NH4_a W.H. sol. NH4 attn. 0.00023 0.00001 0.001 l og 10 u ni WH_Sl_NH4_b W.H. sol. NH4 attn. 0.6 0.3 0.8 u ni WH_Sl_NH4_Cb ppm W.H. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni WH_Sl_ON_a W.H. sol. OrgN attn. 0.00023 0.00001 0.001 l og 10 u ni WH_Sl_ON_b W.H. sol. OrgN attn. 0.6 0.3 0.8 u ni WH_Sl_ON_Cb ppm W.H. sol. OrgN bckgrnd. Con. 1.6 0.5 3 u ni WH_Sl_P_a W.H. sol. P attn. 0.000009 0.000001 0.0001 l og 10 u ni WH_Sl_P_b W.H. sol. P attn. 0.65 0.3 0.8 u ni WH_Sl_P_Cb ppm W.H. sol. P bckgrnd. Con. 0.15 0.005 0.3 l og 10 u ni WH_TSS_a W.H. sol. Sus. Solids attn. 0.0003 0.00003 0.003 l og 10 u ni WH_TSS_b W.H. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni PAGE 179 179 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 WH_TSS_Cb ppm W.H. sol. Sus. Solids bckgrnd. Con. 5 2 15 u ni WH_Sd_NH4_a W.H. sed. NH4 attn. 0.0003 0.00003 0.003 l og 10 u ni WH_Sd_NH4_b W.H. sed. NH4 attn. 0.5 0.3 0.8 u ni WH_Sd_NH4_Cb ppm W.H. sed. NH4 bckgrnd. Con. 0.005 0.001 0.05 l og 10 u ni WH_Sd_ON_a W.H. sed. OrgN attn. 0.0003 0.00003 0.003 l og 10 u ni WH_Sd_ON_b W.H. sed. OrgN attn. 0.5 0.3 0.8 u ni WH_Sd_ON_Cb ppm W.H. sed. OrgN bckgrnd. Con. 0.05 0.015 0.15 l og 10 u ni WH_Sd_P_a W.H. sed. P attn. 0.00003 0.000003 0.0003 l og 10 u ni WH_Sd_P_b W.H. sed. P attn. 0.5 0.3 0.8 u ni WH_Sd_P_Cb ppm W.H. sed. P bckgrnd. Con. 0.05 0.005 0.15 l og 10 u ni WH_GW _ NO3_a W.H. groundwater NO3 attn. 0.8 0.4 0.95 u ni WH_GW _ NH4_a W.H. groundwater NH4 attn. 0.8 0.4 0.95 u ni WH_GW _ ON_a W.H. groundwater OrgN attn. 0.8 0.4 0.95 u ni WH_GW _ P_a W.H. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni CP_Sl_BOD_a Cyp. sol. BOD attn. 0.0004 0.00003 0.003 l og 10 u ni CP_Sl_BOD_b Cyp. sol. BOD attn. 0.5 0.3 0.8 u ni CP_Sl_BOD_Cb ppm Cyp. sol. BOD bckgrnd. Con. 0.5 1 5 u ni CP_Sl_NO3_a Cyp. sol. NO3 attn. 0.00025 0.00001 0.001 l og 10 u ni CP_Sl_NO3_b Cyp. sol. NO3 attn. 0.6 0.3 0.8 u ni CP_Sl_NO3_Cb ppm Cyp. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni CP_Sl_NH4_a Cyp. sol. NH4 attn. 0.00025 0.00001 0.001 l og 10 u ni CP_Sl_NH4_b Cyp. sol. NH4 attn. 0.6 0.3 0.8 u ni CP_Sl_NH4_Cb ppm Cyp. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni CP_Sl_ON_a Cyp. sol. OrgN attn. 0.00025 0.00001 0.001 l og 10 u ni CP_Sl_ON_b Cyp. sol. OrgN attn. 0.6 0.3 0.8 u ni CP_Sl_ON_Cb ppm Cyp. sol. OrgN bckgrnd. Con. 1.6 0.5 3 u ni CP_Sl_P_a Cyp. sol. P attn. 0.000009 0.000001 0.0001 l og 10 u ni CP_Sl_P_b Cyp. sol. P attn. 0.65 0.3 0.8 u ni CP_Sl_P_Cb ppm Cyp. sol. P bckgrnd. Con. 0.15 0.005 0.3 l og 10 u ni CP_TSS_a Cyp. sol. Sus. Solids attn. 0.0004 0.00003 0.003 l og 10 u ni CP_TSS_b Cyp. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni CP_TSS_Cb ppm Cyp. sol. Sus. Solids bckgrnd. Con. 10 2 15 u ni CP_Sd_NH4_a Cyp. sed. NH4 attn. 0.0004 0.00003 0.003 l og 10 u ni CP_Sd_NH4_b Cyp. sed. NH4 attn. 0.5 0.3 0.8 u ni CP_Sd_NH4_Cb ppm Cyp. sed. NH4 bckgrnd. Con. 0.01 0.001 0.05 l og 10 u ni CP_Sd_ON_a Cyp. sed. OrgN attn. 0.0004 0.00003 0.003 l og 10 u ni CP_Sd_ON_b Cyp. sed. OrgN attn. 0.5 0.3 0.8 u ni CP_Sd_ON_Cb ppm Cyp. sed. OrgN bckgrnd. Con. 0.1 0.015 0.15 l og 10 u ni CP_Sd_P_a Cyp. sed. P attn. 0.00004 0.000003 0.0003 l og 10 u ni PAGE 180 180 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 CP_Sd_P_b Cyp. sed. P attn. 0.5 0.3 0.8 u ni CP_Sd_P_Cb ppm Cyp. sed. P bckgrnd. Con. 0.1 0.005 0.15 l og 10 u ni CP_GW _ NO3_a Cyp. groundwater NO3 attn. 0.8 0.4 0.95 u ni CP_GW _ NH4_a Cyp. groundwater NH4 attn. 0.8 0.4 0.95 u ni CP_GW _ ON_a Cyp. groundwater OrgN attn. 0.8 0.4 0.95 u ni CP_GW _ P_a Cyp. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni WF_Sl_BOD_a W.M.F. sol. BOD attn. 0.0004 0.00003 0.003 l og 10 u ni WF_Sl_BOD_b W.M.F. sol. BOD attn. 0.5 0.3 0.8 u ni WF_Sl_BOD_Cb ppm W.M.F. sol. BOD bckgrnd. Con. 0.5 1 5 u ni WF_Sl_NO3_a W.M.F. sol. NO3 attn. 0.00023 0.00001 0.001 l og 10 u ni WF_Sl_NO3_b W.M.F. sol. NO3 attn. 0.6 0.3 0.8 u ni WF_Sl_NO3_Cb ppm W.M.F. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni WF_Sl_NH4_a W.M.F. sol. NH4 attn. 0.00023 0.00001 0.001 l og 10 u ni WF_Sl_NH4_b W.M.F. sol. NH4 attn. 0.6 0.3 0.8 u ni WF_Sl_NH4_Cb ppm W.M.F. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni WF_Sl_ON_a W.M.F. sol. OrgN attn. 0.00023 0.00001 0.001 l og 10 u ni WF_Sl_ON_b W.M.F. sol. OrgN attn. 0.6 0.3 0.8 u ni WF_Sl_ON_Cb ppm W.M.F. sol. OrgN bckgrnd. Con. 1.5 0.5 3 u ni WF_Sl_P_a W.M.F. sol. P attn. 0.000009 0.000001 0.0001 l og 10 u ni WF_Sl_P_b W.M.F. sol. P attn. 0.65 0.3 0.8 u ni WF_Sl_P_Cb ppm W.M.F. sol. P bckgrnd. Con. 0.15 0.005 0.3 l og 10 u ni WF_TSS_a W.M.F. sol. Sus. Solids attn. 0.0004 0.00003 0.003 l og 10 u ni WF_TSS_b W.M.F. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni WF_TSS_Cb ppm W.M.F. sol. Sus. Solids bckgrnd. Con. 5 2 15 u ni WF_Sd_NH4_a W.M.F. sed. NH4 attn. 0.0004 0.00003 0.003 l og 10 u ni WF_Sd_NH4_b W.M.F. sed. NH4 attn. 0.5 0.3 0.8 u ni WF_Sd_NH4_Cb ppm W.M.F. sed. NH4 bckgrnd. Con. 0.005 0.001 0.05 l og 10 u ni WF_Sd_ON_a W.M.F. sed. OrgN attn. 0.0004 0.00003 0.003 l og 10 u ni WF_Sd_ON_b W.M.F. sed. OrgN attn. 0.5 0.3 0.8 u ni WF_Sd_ON_Cb ppm W.M.F. sed. OrgN bckgrnd. Con. 0.05 0.015 0.15 l og 10 u ni WF_Sd_P_a W.M.F. sed. P attn. 0.00004 0.000003 0.0003 l og 10 u ni WF_Sd_P_b W.M.F. sed. P attn. 0.5 0.3 0.8 u ni WF_Sd_P_Cb ppm W.M.F. sed. P bckgrnd. Con. 0.05 0.005 0.15 l og 10 u ni WF_GW _ NO3_a W.M.F. groundwater NO3 attn. 0.8 0.4 0.95 u ni WF_GW _ NH4_a W.M.F. groundwater NH4 attn. 0.8 0.4 0.95 u ni WF_GW _ ON_a W.M.F. groundwater OrgN attn. 0.8 0.4 0.95 u ni WF_GW _ P_a W.M.F. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni FM_Sl_BOD_a F.M. sol. BOD attn. 0.0005 0.00003 0.003 l og 10 u ni FM_Sl_BOD_b F.M. sol. BOD attn. 0.5 0.3 0.8 u ni PAGE 181 181 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 FM_Sl_BOD_Cb ppm F.M. sol. BOD bckgrnd. Con. 0.5 1 5 u ni FM_Sl_NO3_a F.M. sol. NO3 attn. 0.0003 0.00001 0.001 l og 10 u ni FM_Sl_NO3_b F.M. sol. NO3 attn. 0.6 0.3 0.8 u ni FM_Sl_NO3_Cb ppm F.M. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni FM_Sl_NH4_a F.M. sol. NH4 attn. 0.0003 0.00001 0.001 l og 10 u ni FM_Sl_NH4_b F.M. sol. NH4 attn. 0.6 0.3 0.8 u ni FM_Sl_NH4_Cb ppm F.M. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni FM_Sl_ON_a F.M. sol. OrgN attn. 0.0003 0.00001 0.001 l og 10 u ni FM_Sl_ON_b F.M. sol. OrgN attn. 0.6 0.3 0.8 u ni FM_Sl_ON_Cb ppm F.M. sol. OrgN bckgrnd. Con. 1.9 0.5 3 u ni FM_Sl_P_a F.M. sol. P attn. 0.00001 0.000001 0.0001 l og 10 u ni FM_Sl_P_b F.M. sol. P attn. 0.65 0.3 0.8 u ni FM_Sl_P_Cb ppm F.M. sol. P bckgrnd. Con. 0.05 0.005 0.3 l og 10 u ni FM_TSS_a F.M. sol. Sus. Solids attn. 0.0005 0.00003 0.003 l og 10 u ni FM_TSS_b F.M. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni FM_TSS_Cb ppm F.M. sol. Sus. Solids bckgrnd. Con. 5 2 15 u ni FM_Sd_NH4_a F.M. sed. NH4 attn. 0.0005 0.00003 0.003 l og 10 u ni FM_Sd_NH4_b F.M. sed. NH4 attn. 0.5 0.3 0.8 u ni FM_Sd_NH4_Cb ppm F.M. sed. NH4 bckgrnd. Con. 0.005 0.001 0.05 l og 10 u ni FM_Sd_ON_a F.M. sed. OrgN attn. 0.0005 0.00003 0.003 l og 10 u ni FM_Sd_ON_b F.M. sed. OrgN attn. 0.5 0.3 0.8 u ni FM_Sd_ON_Cb ppm F.M. sed. OrgN bckgrnd. Con. 0.05 0.015 0.15 l og 10 u ni FM_Sd_P_a F.M. sed. P attn. 0.00005 0.000003 0.0003 l og 10 u ni FM_Sd_P_b F.M. sed. P attn. 0.5 0.3 0.8 u ni FM_Sd_P_Cb ppm F.M. sed. P bckgrnd. Con. 0.05 0.005 0.15 l og 10 u ni FM_GW _ NO3_a F.M. groundwater NO3 attn. 0.8 0.4 0.95 u ni FM_GW _ NH4_a F.M. groundwater NH4 attn. 0.8 0.4 0.95 u ni FM_GW _ ON_a F.M. groundwater OrgN attn. 0.8 0.4 0.95 u ni FM_GW _ P_a F.M. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni WP_S l _ BOD_a W.P. sol. BOD attn. 0.0004 0.00003 0.003 l og 10 u ni WP_S l _ BOD_b W.P. sol. BOD attn. 0.5 0.3 0.8 u ni WP_S l _ BOD_Cb ppm W.P. sol. BOD bckgrnd. Con. 0.5 1 5 u ni WP_S l _ NO3_a W.P. sol. NO3 attn. 0.00012 0.00001 0.001 l og 10 u ni WP_S l _ NO3_b W.P. sol. NO3 attn. 0.6 0.3 0.8 u ni WP_S l _ NO3_Cb ppm W.P. sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni WP_S l _ NH4_a W.P. sol. NH4 attn. 0.00012 0.00001 0.001 l og 10 u ni WP_S l _ NH4_b W.P. sol. NH4 attn. 0.6 0.3 0.8 u ni WP_S l _ NH4_Cb ppm W.P. sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni WP_S l _ ON_a W.P. sol. OrgN attn. 0.00012 0.00001 0.001 l og 10 u ni WP_S l _ ON_b W.P. sol. OrgN attn. 0.6 0.3 0.8 u ni PAGE 182 182 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 WP_S l _ ON_Cb ppm W.P. sol. OrgN bckgrnd. Con. 1.6 0.5 3 u ni WP_S l _ P_a W.P. sol. P attn. 0.000004 0.000001 0.0001 l og 10 u ni WP_S l _ P_b W.P. sol. P attn. 0.65 0.3 0.8 u ni WP_S l _ P_Cb ppm W.P. sol. P bckgrnd. Con. 0.15 0.005 0.3 l og 10 u ni WP_TSS_a W.P. sol. Sus. Solids attn. 0.0004 0.00003 0.003 l og 10 u ni WP_TSS_b W.P. sol. Sus. Solids attn. 0.5 0.3 0.8 u ni WP_TSS_Cb ppm W.P. sol. Sus. Solids bckgrnd. Con. 5 2 15 u ni WP_S d _ NH4_a W.P. sed. NH4 attn. 0.0004 0.00003 0.003 l og 10 u ni WP_S d _ NH4_b W.P. sed. NH4 attn. 0.5 0.3 0.8 u ni WP_S d _ NH4_Cb ppm W.P. sed. NH4 bckgrnd. Con. 0.005 0.001 0.05 l og 10 u ni WP_S d _ ON_a W.P. sed. OrgN attn. 0.0004 0.00003 0.003 l og 10 u ni WP_S d _ ON_b W.P. sed. OrgN attn. 0.5 0.3 0.8 u ni WP_S d _ ON_Cb ppm W.P. sed. OrgN bckgrnd. Con. 0.05 0.015 0.15 l og 10 u ni WP_S d _ P_a W.P. sed. P attn. 0.00004 0.000003 0.0003 l og 10 u ni WP_S d _ P_b W.P. sed. P attn. 0.5 0.3 0.8 u ni WP_S d _ P_Cb ppm W.P. sed. P bckgrnd. Con. 0.05 0.005 0.15 l og 10 u ni WP_GW _ NO3_a W.P. groundwater NO3 attn. 0.8 0.4 0.95 u ni WP_GW _ NH4_a W.P. groundwater NH4 attn. 0.8 0.4 0.95 u ni WP_GW _ ON_a W.P. groundwater OrgN attn. 0.8 0.4 0.95 u ni WP_GW _ P_a W.P. groundwater P attn. 0.04 0.01 0.1 l og 10 u ni D_Sl_BOD_a Dryland sol. BOD attn. 0.0001 0.00001 0.001 l og 10 u ni D_Sl_BOD_b Dryland sol. BOD attn. 0.7 0.45 0.85 u ni D_Sl_BOD_Cb ppm Dryland sol. BOD bckgrnd. Con. 0.5 0.1 5 l og 10 u ni D_Sl_NO3_a Dryland sol. NO3 attn. 0.0001 0.00001 0.001 l og 10 u ni D_Sl_NO3_b Dryland sol. NO3 attn. 0.7 0.45 0.85 u ni D_Sl_NO3_Cb ppm Dryland sol. NO3 bckgrnd. Con. 0.02 0.001 0.1 l og 10 u ni D_Sl_NH4_a Dryland sol. NH4 attn. 0.0001 0.00001 0.001 l og 10 u ni D_Sl_NH4_b Dryland sol. NH4 attn. 0.7 0.45 0.85 u ni D_Sl_NH4_Cb ppm Dryland sol. NH4 bckgrnd. Con. 0.015 0.0015 0.15 l og 10 u ni D_Sl_ON_a Dryland sol. OrgN attn. 0.0001 0.00001 0.001 l og 10 u ni D_Sl_ON_b Dryland sol. OrgN attn. 0.7 0.45 0.85 u ni D_Sl_ON_Cb ppm Dryland sol. OrgN bckgrnd. Con. 0.2 0.05 0.5 l og 10 u ni D_Sl_P_a Dryland sol. P attn. 0.00001 0.000001 0.0001 l og 10 u ni D_Sl_P_b Dryland sol. P attn. 0.6 0.45 0.85 u ni D_Sl_P_Cb ppm Dryland sol. P bckgrnd. Con. 0.1 0.005 0.3 l og 10 u ni D_TSS_a Dryland sol. Sus. Solids attn. 0.0001 0.00001 0.001 l og 10 u ni D_TSS_b Dryland sol. Sus. Solids attn. 0.7 0.45 0.85 u ni D_TSS_Cb ppm Dryland sol. Sus. Solids bckgrnd. Con. 4 2 15 u ni D_Sd_NH4_a Dryland sed. NH4 attn. 0.00001 0.000001 0.0001 l og 10 u ni PAGE 183 183 Table C 1. Continued Parameter Name 1 Units Type 2 Base Value Min Max Distribution 3 D_Sd_NH4_b Dryland sed. NH4 attn. 0.7 0.45 0.85 u ni D_Sd_NH4_Cb ppm Dryland sed. NH4 bckgrnd. Con. 0.007 0.001 0.05 l og 10 u ni D_Sd_ON_a Dryland sed. OrgN attn. 0.00001 0.000001 0.0001 l og 10 u ni D_Sd_ON_b Dryland sed. OrgN attn. 0.7 0.45 0.85 u ni D_Sd_ON_Cb ppm Dryland sed. OrgN bckgrnd. Con. 0.07 0.015 0.15 l og 10 u ni D_Sd_P_a Dryland sed. P attn. 0.00001 0.000001 0.0001 l og 10 u ni D_Sd_P_b Dryland sed. P attn. 0.7 0.45 0.85 u ni D_Sd_P_Cb ppm Dryland sed. P bckgrnd. Con. 0.06 0.005 0.15 l og 10 u ni D_GW _ NO3_a Dryland groundwater NO3 attn. 0.8 0.4 0.95 u ni D_GW _ NH4_a Dryland groundwater NH4 attn. 0.8 0.4 0.95 u ni D_GW _ ON_a Dryland groundwater OrgN attn. 0.8 0.4 0.95 u ni D_GW _ P_a Dryland groundwater P attn. 0.04 0.01 0.1 l og 10 u ni RO_HYD Runoff hydrograph profiles 1 4 disc u ni GW_HYD Groundwater hydrograph profiles 1 4 disc u ni D_RO_vel m/d Dryland runoff velocity 20000 1000 40000 u ni D_RO_kDF d Dryland runoff constant delay 0 0 3 u ni D_GW_vel m/d Dryland groundwater velocity 0 2 200 l og 10 u ni D_GW_kDF d Dryland groundwater constant delay 0 0 15 u ni W_RO_vel m/d Wetland runoff velocity 20000 1000 40000 u ni W_RO_kDF d Wetland runoff constant delay 0 0 3 u ni W_GW_vel m/d Wetland groundwater velocity 0 2 200 l og 10 u ni W_GW_kDF d Wetland groundwater constant delay 0 0 15 u ni dep_err error in depth measurement 20% 20% u ni wid_err error in width measurement 20% 20% u ni n Manning's roughness coefficient 0.05 0.02 0.08 u ni 1 For the water quality attenuation type of parameters first part of the parameter symbol corresponds to land use or reach type, second part indicates if the constituent is in soluble or sediment form, third part indicates the water quality constituent, and the last part indicates the parameter type a, b, and Cb 2 sol. soluble, attn. attenuation, bckgrnd. background, Con. concentration, sed. sediment, sus. suspended, B.S. Bay Swamp, S.L.S. Stream and Lake Swamp, W.H. Wetland Hardwood, Cyp. Cypress, W.M. F. Wetland Mix Forest, F.M. Freshwater Marshes 3 l og10 u ni l oguniform distribution to the base 10, u ni. u niform distribution , disc uni discrete uniform distribution PAGE 184 184 Figure C 1. Runoff unit hydrograph options used for the EE analysis. RO_Hydr ograph 3 is the default runoff hydrogrpah PAGE 185 185 Figure C 2. Sub surface unit hydrograph options used for the EE analysis. Sub surface Unit Hydrograph 1 is the default groundwater unit hydrograph in WAM. Note that all hydrogrpahs are of 90 day duration. For brevity only first 10 days are shown. PAGE 186 186 Figure C RO HYD 1 is the default runoff hydrogrpah PAGE 187 187 LIST OF REFERENCES Allen, R.G., L.S. Pereira, D. Raes, and M. Smith, 1998. Crop Evapotranspirat ion, Guidelines for computing crop water requirements. Irrigation and Drainage Paper 56. Food and Agricultural Organization, Rome. Andres, T.H. and W.C. Hajas, 1993. Using iterated fractional factorial design to screen parameters in sensitivity analysis of a probabilistic risk assessment model. Proceedings of the Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Application, Karlsruhe, Germany. 19 23 April 1993. Arabi, M., R.S. Govindaraju, B. Engel, and M. Hantush, 2007. Multiobjective sensitivity analysis of sediment and nitrogen processes with a watershed model. Water Resources Research 43: W06409. doi:10.1029/2006WR005463 Baginska, B., W . Milne Home, and P. S. Cornish, 2003. Modelling nutrient transport in Currency Creek , NSW, with AnnAGNPS and PEST. Environmental Modelling and Software 18(8 9): 801 808. Bahremand, A. and F. De Smedth, 2008. Distributed hydrologic modeling and sensitivity analysis in a Torysa watershed, Slovakia. Water Resource Manage ment 22: 393 408. Bar oni, G. and S. Tarantola, 2014. A general probabilistic framework for uncertainty and global sensitivity analysis of deterministic models: A hydrologic case study . Environmental Modelling and Software 51: 26 34. Bedford, T. , 1998. Sensitivity indices for (Tree) dependent variables. In: Proceedings of the second international symposium on sensitivity analysis of model output. Venice, Italy. Benke, K. K., K. E. Lowell, and A. J. Hamilton, 2008. Parameter uncertainty, sensitivity analysis and prediction error in a water balance hydrological model. Math ematical and Computer Modell ing 47(11 12): 1134 1149. Bennett, N.D., B.F.W. Croke, G. Guariso, J.H.A. Guillaume, S.H. Hamilton, A.J. Jakeman, S. Marsili Libelli, L.T.H. Ne wham, J.P. Norton, C. Perrin, S.A. Pierce, B. Robson, R. Seppelt, A.A. Voinov, B.D. Fath, and V. Andreassian, 2013. Characterising performance of environmental models. Environmental Mo delling and Software 40: 1 20. d oi.org/10.1016/j.envsoft.2012.09.011 Bergman, M.J., W. Green, and L.J. Donnangelo, 2002. Calibration of storm loads in the South Prong watershed, Florida using BASINS/HSPF. Journal of the AWRA 38(5): 1423 1436. Bettonvil, B. and J.P.C. Kleijnen, 1997. Searching for important factors in simulation models with many factors: sequential bifurcation. European Journal of Operation Research 96: 180 194. doi.org/10.1016/S0377 2217(96)00156 7 Bettonvil, B. , 1990. Detection of important factors by sequenti al bifurcation, Tilburg University Press, Tilburg. PAGE 188 188 Beven , K.J., P.J. Smith, and A. Wood, 2011. On the colour and spin of epistemic error (and wha t we might do about it). Hydrological Earth System Science Discussion 8: 5355 5386. d oi:10.51194/hessd 8 5355 2 011 Beven, K. , 2006. On undermining the science? Hydrological Processes 20(14): 3141 3146. Beven, K. , 2001. Rainfall runoff modeling The Primer, John Wiley and Sons, Chichester, England, ISBN 978 0 471 98552 2. Beven, K. and A. Binley, 1992. The future of distributed models: model calibration and uncertainty pre diction. Hydrological Processes 6: 279 298. Beven, K. , 1989. Changing ideas in hydrology: The case of physically based models. Journal of Hydrology 105(1 2): 157 172. Bicknell, B .R., J. C . Imhoff, J. L. Kittle, Jr., T.H. Jobes, and A. S. Donigian, Jr., 2001. Hydrological Simulation Program Fortran, Version 12: User's Manual. Mountain View, Cal.: Aqua Terra Consultants. Bingner, R.L. and F. D. Theurer, 2001. AnnAGNPS Technical Proces ses: Documentation Version 2. Available at: www.sedlab.olemiss.edu/AGNPS.html . Blasone, R. S., J.A. Vrugt, H. Madsen, D. Rosbjerg, B.A. Robinson, and G. A. Zyvoloski , 2008. Generalized likelihood unce rtainty estimation (GLUE) using adaptive Markov Chain Monte Carlo sampling. Advances in Water Resources 31: 630 648. Borah, D.K. and M. Bera , 2004. Watershed scale hydrologic and nonpoint source pollution models: Review of applications. Transactions of ASABE 47(3): 789 803 . Borah, D. K. and M. Bera , 2003. Watershed scale hydrologic and nonpoint source pollution models: Review of mathematical bases. Transactions of ASAE 46(6): 1553 1566. Bottcher, A.B., B.J. Whiteley, James, A.I., and J.G. Hiscock, 2012. Watershed Assessment Model (WAM): Model use, calibration and v alidation. Transactions of ASABE 55(4): 1367 1383. Bottcher, A.B., J.G. Hiscock, a nd B.M. Jacobson , 2002. WAM View, A GIS Approach to Watershed Assessment Modeling. Proceeding of the Watershed 2 002 Pre Conference Modeling Workshop. Ft. Lauderdale, FL. Publisher: Water Environment Federation. Alexandria, VA . Bottcher, A.B., N.B. Pickering, and A. B. Cooper. 1998. EAAMOD FIELD: A flow and phosphorus model for high water tables. In Proc. 7th Annual D rainage Symposium, 599 606. St. Joseph, Mich.: ASAE. Bouraoui, F., I. Braud, and T . A. Dillaha, pollution model for water, sediment, and nutrient losses. In Mathematical Models of PAGE 189 189 Small Watershed Hydrology and Ap eds. Highlands Ranch, Colo.: Water Resources Publications . Bristow, K.L. and G.S. Campbell, 1984. On the relationship between incoming solar radiation and daily maximum and minimum temperature. Agricultural Forest Meteorology 31: 159 166. Brown, M.T., 1995. Water quality and wildlife values of wetlands. Center for Wetlands, University of Florida , Gainesville, FL. Byrne, M.J. and M.S. Wood, 2012. Concentrations and loads of nutrients in the tributaries of the Lake Okeechobee Watershed, South Central Florida, Water Years 2004 2008. USGS Reston, VA. Cacuci, D.G. and M. Ionescu Bujor, 2004. A compara tive review of sensitivity and uncertainty analysis of large scale systems II: Statistical Methods. Nuclear Science and Engineering 147(3): 204 217. Campolongo, F., A. Saltelli, and J. Cariboni, 2011. From screening to qualitative sensitivity analysis. A u nified approach. Computer Physics Communications 182: 978 988. doi:10.1016/j.cpc.2010.12.039 Campolongo, F., J. Cariboni, and A. Saltelli, 2007. An effective screening design for sensitivity analysis of large models. Environmental Modelling and Software 22 : 1509 1518. doi:10.1016/j.envsoft.2006.10.004 Campolongo, F., S. Tarantola , and A. Saltelli, 1999. Tackling quantitatively large dimensionality problems. Computer Physics and Communications 117: 75 85. doi.org/10.1016/S0010 4655(98)00165 9 . Campolongo, F. and A. Saltelli, 1997. Sensitivity analysis of an environmental model: an application of different analysis methods. Reliability Engineering and System Safety 57(1): 49 69. doi.org/10.1016/S0951 8320(97)00021 5 Cariboni, J., D. Gatelli, R. Liska, and A. S altelli, 2007. The role of sensitivity analysis in ecological modelling. Ecological Modelling 203: 167 182. doi:10.1016/j.ecolmodel.2005.10.045. Carrow, R.N. , 1995. Drought resistance aspect of turfgrasses in the southeast: Evapotranspiration and crop coef ficients. Crop Science 35: 1685 1690. Chaubey, I., A.S. Cotter, T.A. Costello, and T. S. Soerens, 2005. Effect of DEM data resolution on SWAT output uncertainty. Hydrol ogical Processes 19: 621 628. Chaubey, I., C.T. Haan, S. Grunwald, and J. M. Salsbury, 199 9. Uncertainty in model parameters due to spatial variability of rainfall. Journal of Hydrology. 220: 48 61. PAGE 190 190 Chebud, Y., G.M. Naja, and R. Rivero, 2011. Phosphorus runoff assessment in a watershed. Journal of Environmental Monitoring 13: 66 73. Cho, J., D. Bosch, R. Lowrance, T. Strickland, and G. Vellidis, 2009. Effect of spatial distribution of rainfall on temporal and spatial uncertainty o f SWAT output. Transactions of ASABE 52(5): 1545 1555. Chow, V.T., 1959 . Open channel hydraulics: New York, McGraw Hi ll, 680 p. Christiaens, K., 2002. Use of sensitivity and uncertainty measures in distributed hydrological modeling with an application to the MIKE SHE model Water Resources Research 38(9): 1169 . doi:10.1029/2001WR000478 Cibin, R., K.P. Sudheer, and I. Cha ubey , 2010. Sensitivity and identifiability of stream flow gene ration parameters of SWAT model. Hydrological Processes 24(9): 1133 1148. doi 10.1002/hyp.7568. Ciric, C., P. Ciffroy, and S. Charles, 2012. Use of sensitivity analysis to identify influential and non influential parameters within an aquatic ecosystem model. Ecological Modelling 246: 119 130. doi.org/10.1016/j.ecolmodel.2012.06.024 Cloke, H., F. Pappenberger, and J. P. Renaud, 2008. Multi method global sensitivity analysis (MMGSA) for modelling floodplain hydrological processes. Hydrological Processes 22: 1660 1674. Confalonieri, R., G. Bellocchi, S. Bregagl io, M. Donatelli, and M. Acutis, 2010. Comparision of sensitivity techniques: A case study with the rice model WARM. Ecological Modeling 221: 1897 1906. Cotter, A.S., I. Chaubey, T. A. Costello, T.S. Soerens, and M. A. Nelson, 2003. Water quality model output uncertainty as affected by spat ial resolution of input data. Journal of the AWRA 39(4): 977 986. Crosetto, M. and S. Tarantola , 2001. Uncertainty and sensitivity analysis: tools for GIS based model implementation. International Journal of Geop hysical Information and Science 15(5): 415 437. Cukier, R.I., J.H. Schaibly, and K.E. Schuler, 1978. Nonlinear sensitivity analysis of multi parameter model systems. Journal of Computational Physics 26 (1):1 42. Cukier, R.I., C.M. Fortuin, K.E. Schuler, A.G. Petschek, and J.H. Schaibly, 1973. Study of sensitivity of coupled reaction systems to uncertainties in rate coeffici ents. I Theory. Journ al of Chemical Physics 59: 3873 3878. Cullmann, J., V. Mishra, and R. Peters, 2006. Flow analysis with WaSiM ETH model parameter sensi tivity at different scales Advances in Geosci ences 9: 73 77. PAGE 191 191 Daniel, E.B., J.V. Camp, E.J. LeBoeuf, J.R. Penrod, J. P. Dobb ins , and M. D. Abkowitz , 2011. Watershed modeling and its applicat ions: A state of the art review. The Open Hydrology Journal 5: 26 50 . Dardanelli, J.L., J.T. Ritchie, M. Calmon, J.M. Andriani, and D.J. Collino, 2004. An empirical model for root water uptak e. Field Crop Research 87: 59 71. Dardanelli, J.L., O.A. B achmeier, R. Sereno, and R. Gil, 1997. Rooting depth and soil water extraction patterns of different crops in a silty loam Haplustoll. Field Crops Research 54: 29 38. Debarry, P.A. and R.G. Quimpo, 1999. GIS modules and distributed models of the watershed. ASCE. ISBN: 0 7844 0443 7. Deflandre, A., R. J. Williams , F.J. Elorza, J. Mira, and D. B. Boorman, 2006. Analysis of the QUESTOR water quality model using a Fourier amplitude sensitivity test (FAST) for two UK rivers . Science of the Total Environment 360(1 3): 290 304. Demaria, E. M., B. Nijssen, and T. Wagener, 2007. Monte Carlo sensitivity analysis of land surface parameters using the variable i nfiltration capacity model. Journal of Ge ophys ical Research 112: D 11113, doi:10.1029/2006JD007534 Dimov, I. and R. Georgieva, indices. Mathematics and Computers in Simulation 81: 506 514. European Comission (EC) , 2009. Impact Assessment Guidelines. Technical Repor t 92. http://ec.europa.eu/smartregulation/impact/commission_guidelines/docs/iag_2009_en.pdf (accessed: November 2013) Fang, K. T., D.K. Lin, P. Winker, and Y. Zhang, 2000. Uniform design: Theory and application. Technometrics 42(3): 237 248. Florida Department of Environemental Protection, 2005. Nutrient Total Maximum Daily Loads for the Lake Lena. October 2005. Florida Department of Transportation, 2004. Drainage Hankbook Hydrology. Tallahassee, FL. Foglia, L., M.C. Hill, S. W. Mehl, and P. Burlando, 2009. Sensitivity analysis, calibration, and testing of a distributed hydrological model using error based weighting and one o bjective fu nction, Water Resources Research 45: W 06427. doi:10.1029/2008WR007255 Fraisse, C.W., E.M. Gelcer, and P. Woli, 2011. Drought decision support tool: Introducing the Agricultural Reference Index for Drought ARID, University of Florida IFAS Extension, Pub lication # AE 469. http://edis.ifas.ufl.edu/ae469 . (a ccessed on 11 April, 2012) Francos, A., F.J. Elorza, F. Bouraou i, G. Bidoglio, and L. Galbiati, 2003. Sensitivity analysis of distributed environmental simulation models: understanding the model behavior in PAGE 192 192 hydrological studies at the catchment scale. Reliability Engineering and Sys tems Safety 79: 205 218. Frere, M. and G. Popov, 1986. Early agromet eorological crop yield assessment. Plant Production and Protection Paper No. 73. Rome, Italy: Food and Agricultural Organization of the United Nations. Gallagher, M. and J. Doherty , 2007. Parameter estimation and uncertainty analysis for a watershed model. Environ mental Modell ing and Soft ware 22: 1000 1020. Gan, Y., Q. Duan, W. Gong, C. Tong, Y. Sun, W. Chu, A. Ye, C. Miao, and Z. Di, 2014. A comprehensive evaluation of various sensitivity analysis methods: A case study with a hydrologic model. Environmenta l Modelling and Software 51: 269 285. Ge, Q., B. Ciuffo, and M. Menendez, 2014. Combining screening and metamodel based methods: an efficient sequential approach for the sensitivity analysi s of model outputs. Reliability Engineering and System Safety (in p ress). Gijsman, A. J., S.S. Jagtap, and J.W. Jones, 2003. Wading through a swamp of complete confusion: how to choose a method for estimating soil water retention parameters for crop models. European Journal of Agronomy 18: 77 106. Govindaraju, R.S. and A.R . Rao, 2000. Artificial Neural Networks in Hydrology. Kluwer Academic Publishers, Boston , MA . Graham, W.D., A.S. Carpena, W. Skaggs, and A. Shirmohammadi, 2009. Peer review of the Watershed Assessment Model, Final Report. Greets, S., D. Raes, M. Garcia, R. Miranda, J.A. Cusicanqui, C. Taboada, J. Mendoza, R. Huanca, A. Mamani, O. Condori, J. Mamani, B. M orales, V. Osco, and P. Steduto, 2009. Simulating yield response of quinoa to water availability with AquaCrop. Agronomy Journal 101(3): 449 508. Gupta, H.V. and H. Kling, 2011. On typical range, sensitivity, and normalization of Mean Squared Error and Nash Sutcliffe Efficiency type metrics. Water Resources Research 47, W10601. http://dx.doi.org/10.1029/2011WR010962. Gupta, H.V., S. Soroos hian, and P.O. Yapo , 1999. Status of automatic calibration for hydrologic models: Comparision with multilevel expert calibration . Journal of Hydrological Engineering 4(2): 135 143 . Guzman, J., A. Shirmohammadi, A. Sadeghi, X. Wang, M.L. Chu, M. Jha, P. Parajuli, Y.P. Khare, and J. Hernandez, 2014. Uncertainty and uncertainty analysis in the context of hydr ologic and water quality models. Trans actions of ASABE (in revision). Haan, P.K. and R. W. Skaggs , 2003. Effect of parameter uncertainty on DRAINMOD pr edictions : I. Hydrology and yield. Transactions of ASAE 46(4): 1061 1067. PAGE 193 193 Haan, C.T., D.E. Storm, T. Al Issa, S. Prabhu, G.J. Sabbagh, and D.R. Edwards, 1998. Effect of parameter distribution on uncertainty analy sis of hydrologic models. Transactions of AS AE 41(1): 65 70. Haan , C.T., B. Allred, D.E. Strom, G.J. Sabbagh, and S. Prabhu, 1995. Statisitcal procedure for evaluating Hydrolog ic/ Water Quality models. Transactions of ASAE 38(3): 725 733. Hamby, D.M. , 1994. A review of techniques for parametric sensitivity analysis of environmental models. Environmental Monitoring and Assessment 32: 135 154. Harivand i A.M., J. Baird, J. Hartin, J. M. Henry, and D. Shaw , 2009. Managing turfgrass during drought, Universit y of California, Division of Agricultural and Natural Resources, Publication: 8395, http://anrcatalog.ucdavis.edu/Items/8395.aspx (Accessed on 15 March 2011). Harmel, R.D., P.K. Smith, K.W. Carpena, and B.J. Robson, 2014. Evaluating, interpreting and communicating performance of hydrologic/water quality models considering intended use: A review and recommendations. E nvironmental Modelling and Software 57: 40 51. Harmel, R.D., P.K. Smith, and K.W. Migliaccio, 2009. Modifying goodness of fit indicators to incorporate both measurement and model uncertainty in model calibration and validation. Trans actions of ASABE 53(1): 55 63. Harmel, R.D. and P.K. Smith, 2007. Consideration of measurement uncertainty in the evaluation of goodness of fit in hydrologic and water quality modeling. Journal of Hydrology 337: 326 336. Harmel, R.D., R. J. Cooper, R. M. Slade, R.L. Haney, a nd J.G . Arnold, 2006. Cumulative uncertainty i n measured streamflow and water quality data for small waters heds. Transactions of ASABE 49 (3): 689 701. Helsel, D.R. and R.M. Hirsch, 2002. Statistical Methods in Water Resources Techniques of Water Resources Inves tigations, Book 4. United States Geological Survey, 522 p. http://water.usgs.gov/pubs/twri/twri4a3/ (ChapterA3). Helton, J.C. and F.J. Davis , 2002. Illustration of sampling based methods for uncer tai nty and sensitivity analysis. Risk Analysis 22(3): 591 622. Helton, J.C., J.W. Garner, M.G. Mar ietta, R.D. Rechard, D.K. Rudeen, and P.N. Swift, 1993. Uncertainty and sensitivity analysis obtained in a preliminary performance assessment of the waste isolated pilot plant. Nuclear Science and Engineering 114: 286 331. Herman, J.D., J.B. Kollat, P.M. Reed, and T. Wagener, 2013a . From maps to movies: high resolution time varying sensitivity analysis for spatially distributed watershed model. Hydrology Earth System Science 17: 5109 5126. doi: 10.5194/hess 17 5109 2013 PAGE 194 194 Herman, J.D., J.B. Kollat, P.M. Reed, and T. Wagener, 2013b. Technical note: Method of Morris effectively reduces the computational demands of global sensitivity analysis for distributed watershed models. Hydrology and Earth System Sciences 17: 2893 2903. doi:10.5194/hess 17 2893 2013 Hornberger, G.M. and R.C. Spear, 1981. An approach to th e preliminary analysis of environmental systems. Journal of Environmental Management 12: 7 18. Hosseinipour, E., 2006. Integrated hydrologic modeling for sustainable development and management of urban water supplies, In Urban Groundwater Management and Su stainability, Ed. J.H. Tellam. Springer. High Performance Computing Center, University of Florida, Gainesville, FL. http://www.hpc.ufl.edu/ (accessed in March 2013) Ibrekk, H. and M.G. Morgan, 1987. Graphical comm unication of uncertain quantities to non technical people. Risk Analysis 7: 519 529. Idso, S. B., R.D. Jackson, P.J. Pinter Jr., R .J. Reginato, and J.L. Hatfield, 1981. Normalizing the stress degree day parameter for environmental variability. Agricultural Meteorology 24: 45 55. Iman, R.L., M .E. Johnson, and T.A. Schroeder, 2002. Assessing hurricane effects Part I: Sensitivity analysis. Reliability Engineering and System Safety 78: 131 145. Iman, R.L. and J.C. Helton , 1988. An investigation of uncertainty a nd sensitivity analysis for computer models. Risk Analysis 8(1): 71 90. Iman, R.L., M.J. Shortencarier, and J.D. Jonson, calculation of partial correlation function and standardized coefficients, In: NUREG/CR 4122, SAND85 0044. Albuquerque, NM: Sandia National Laboratories. Iman, R.L. and W.J. Conover , 1982. A distribution free approac h to inducing rank correlation among input variables. Communications in Statistics Simulation and Computati on 11 (3): 311 334. Jaber, F.H. and S. Shukla, 2007. Accuracy of hydrodynamic modeling of flood deterntion reservoirs. Journal of Hydrologic Engineer ing 12(2): 225 230. Jacobson, B. and A. B. Bottcher , 1998. Unique routing algorithm for watershed assessment modeling. ASAE Paper No. 982237. St. Joseph, Mich.: ASAE. Jain, S.K. and K.P. Sudheer, 2008. Fitting of hydrologic models: a close look at the Nash Sutcliffe Index. Journal of Hydrologic Engineering 13: 981 986. Jakeman, A.J., R.A. Letcher, and J.P. Norton, 2006. Ten iterative steps in dev elopment of and evaluation of environmental models. Environmental Modelling and Software 21: 602 614. doi:10.1016/j.envsoft.2006.01.004 . PAGE 195 195 Jansen, M. , 1999. Analysis of variance designs for model output. Computer Physics Communications 117: 35 43. Jawitz, J.W., R. MuÃ±oz Carpena, S. Muller, K.A. Grace, K.A., and A.I. James , 2008. Development, Testing, and Sensitivity an d Uncertainty Analyses of a Transport and Reaction Simulation Engine (TaRSE) for Spatially Distributed Modeling of Phosphorus in South Florida Peat Marsh Wetlands. http://pubs.usgs. gov/sir/2008/5029/pdf/sir2008 5029.pdf (accessed July 2013) Jin, X., C . Xu, Q. Z hang, and V. P. Singh, 2010. Parameter and modeling uncertainty simulated by GLUE and a formal Bayesian method for a conceptual hydrologic model. Journal of Hydrology 383: 147 155. Johnson, M.E., L.M. Moore, and D. Ylvisaker, 1990. Minimax and maximin distance designs. Journal of Statistical Planning and Inference 26: 131 148. Joint Research Center, European Commission, (JRC, EC) Ispra, VA, Italy. http://ipsc.jrc.ec.europa.eu/?id=756 (accessed in December 2012). Jones, J.W., G. Hoogenboom, C.H. Porter, K.J. Boote, W.D. Batchelor, L.A. Hunt, P.W. Wilkens, U. Singh, A.J. Gijsman, and J.T. Ritchie, 2003. The DSSAT croppi ng system model. European Journal of Agronomy 18: 235 265. Kadlec, R.H. and S.D . Wallace , 2009. Treatment Wetlands , Second Edition. CRC Press. Kadlec, R.H. and R.L. Knight, 1996. Treatment Wetlands, First Edition. CRC Press. Kang, W. J. and D. Gilbert, 201 3. Nutrient and Dissolved Oxygen TMDLs for the Lake Jackson. Prepared for Florida Department of Environmental Protection, Tallahassee, FL. Keener, V. W., K.T. Ingram, B. Jacobson, and J.W. Jones, 2007. Oscillation on simulated Ph osphorus loading in South Florida. Transactions of ASABE 50(6) 2081 2089. Kisekka, I., Carpena, Y.P. Khare, and T.H. Boyer, 2013. Sensitivity analysis and parameter estimation for an approximate analytical model of canal aquifer interaction applied in t he C 111 basin. Transactions of ASABE 56(3): 977 992. doi: 10.13031/trans.56.10037 Knisel, W.G. (Ed). 1993. GLEAMS: Groundwater Loading Effects of Agricultural Management Systems. Univ. of Georgia, Coastal Plain Expt. Sta, Bio. and Agri. Engr. Dept., Publ. No. 5 . Kruskal, W.H. and W.A. Wallis, 1952. Use of ranks in on e criterion variance analysis. Journal of the American Statistical Association 47 (260): 583 621. doi:10.1080/01621459.1952.10483441 PAGE 196 196 Kucherenko, S., S. Tarantola, and P . Annoni, 2012. Estimation of global sensitivity indices for models with dependent variables. Computer Physics Communications 183: 937 946. Kucherenko, S., M. Rodriguez Fernan dez, C. Pantelides, and N. Shah, 2009. Monte Carlo evaluation of derivative based global sensitivity measures. Reliability Engineering and System Safety 94: 1135 1148. Legates, D.R. of hydrologic and hydroclimatic model validation. Water Resources Research 35: 23 3 241. Liu, Y., H. Gupta, E. Spr inger, and T. Wagener, 2008. Linking science with environmental decision making: Experiences from an integrated modelling approach to support sustainable water resources management. Environmental Modelling and Software 23: 846 858. doi: 10.1016/j.envsoft.2007.10.007 Loo svelt, L., V.R.N. Pauwels, W.M. Cornelis, G.J.M. De Lannoy, and N.E.C. Verhoest, 2011. Impact of soil hydraulic parameter uncertainty on soil moisture modeli ng. Water Resources Research 47: W03505, doi:10.102 9/2010WR009204. Makler Pick, V., G. Gal, M. Gorfin, M.R. Hipsey , and Y. Carmel, 2011. Sensitivity analysis for complex ecological models A new approach. Environmental Modelling and Software 26: 124 134. doi: 10.1016/j.envsoft.2010.06.010. Manache, G. and C.S. Melching, 2008. Identification of reliable regression and correlation based sensitivity measures for importance ranking of water quality model parameters . En vironmental Mode l ling and Software 23: 549 562. Massey, F. J. , 1951. The Kolmogorov Smirnov test for goodness of fit. Journal of American Statistical Association 46(253): 68 78. McCormick, P., R. T. James, and J. Zhang, 2010. Chapter 10: Lake Okeechobee protection program th Florida Environmental Report. McCuen, R.H., Z. Knight, and A.G. Cutter, 2006. Evaluation of the Nash Sutcliffe Efficiency Index. J ournal of Hydrologic Engineering 11: 597 602. McIntyre, N.R., T. Wagener, H.S. Wheater, and S. C. Chapra, 2003. Risk based m odelling of surface water quality: A case study of t he Charles River, Massachusetts. Journal of Hydrology 274(1 4): 225 247. McKay, M.D., R.J. Beckman, and W.J. Conover , 1979. A comparison of three methods of selecting values of input variables in the anal ysis of output from a computer code. Technomatrics 21: 239 245. McKee, T. B., N.J. Doesken, and J. Kleist, 1993. The relationship of drought frequency and duration to time scales. p. 179 184. In Proc. Conf. on Applied Climatology, 8 th , Anaheim, CA. 17 22 Ja n. 1993. American Meteorological Society, Boston , MA . PAGE 197 197 Meink e, H., G.L. Hammer, and P. Want, 1993. Potential soil water extraction by sunflower on a range of soils. Field Crops Research 32:59 81. Melching, C.S. , 1995. Reliability estimation, Chapter 3 in Computer Models of Watershed Hydrology, ed. V.P. Singh, Water Resources Publications, Littlenton, CO. Menne, M.J., C. N. Williams, Jr., and R.S. Vose, 2012. United States Historical Climatology Network Daily Tempe rature, Precipitation, and Snow Data. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, Oak Ridge, T N . Mi gliaccio, K. W. and I. Chaubey, 2008. Spatial distributions and stochastic parameter influence on SWAT flow and sediment predic tion. Journal of Hydrol ogic Engineering 13(4): 258 269. Migliaccio, K.W. and P. Srivastava. 2007. Hydrologic components of watershe d sca le models. Transactions of ASABE 50(5): 1695 1703. Miller, S.N., D.P. Guertin, and D. C. Go odrich, 2007. Hydrologic model ing uncertainty resulting from land cover misclassification. J ournal of the AWRA 43(4): 1065 1075. Mirchi, A., D. Watkins, and K. Madani , 2009. Chapter 6: Modeling for watershed planning, management and decision making. In Watersheds: Management, restorati on and environmental impact, Editor: J.C. Vaughn, ISBN: 978 1 61668 667 3, Nova Science Publishers Inc. Mishra, A. , 2011. Estimating uncertainty in HSPF based Water Quality Model: Application of Monte Carlo based techniques, Ph.D. Dissertation, Dept. of Bi ological Systems Engineering, Virginia Polytechnic Institute and State University , VA . Mishra, A.K. and V.P . Singh, 2010. A review of drought concepts. Journal of Hydrology 391: 202 216. Moreau, P., V. Viaud, V. Parnaudeau, J. Salmon Monviola, and P. Durand, 2013. An approach for global sensitivity analysis of a complex environmental model to spatial inputs and parameters: A case study of an agro hydrological model. Environmental Modelling and Software 47: 74 87. doi.org/10.1016/j.envsoft.2013.04.00 6. Moriasi, D.N., J.G. Arnold, M.W. Van Liew, R.L. Bingner, R.D. Harmel, and T.L. Veith, 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50(3): 885 900. Morris, M. D. , 1991. Factorial sampli ng plans for preliminary computational experiments. Technometric, 33(2): 161 174. Muleta, M.K., J. W. Nicklow, and E. G. Bekele, 2007. Sensitivity of a distributed watershed simulation model to spatial scale. Journal of Hydrological Engineering 12(2): 163 1 72. PAGE 198 198 Muleta, M. and J. W. Nicklow, 2005. Sensitivity and uncertainty analysis coupled with automatic calibration for a distributed watersh ed model. Journal of Hydrology 306(1 4): 127 145. doi:10.1016/j.jhydrol.2004.09.005. Carp ena, R., Z. Zajac, and Y.M. Kuo, 2007. Global sensitivity and uncertainty analyses of the water quality model VFSMOD W. Transactions of ASABE 50(5): 1719 1732. Carpena, R., G. Vellidis, A. Shirmohammadi , and W.W. Wallender , 2006. Evaluation of modeling tools for TMDL deve lopment and implementation. Transactions of ASABE 49(4): 961 965. Murray, J., B. Brooks Solveson, J. da Forta, and J. Key, 2004. Scientific evaluation and simulation of the six mile Cypress Watershed. Prepared for Lee Scounty, FL and South Florida Water Ma nagement District (SFWMD), FL. Nar asimhan, B. and R. Srinivasan, 2005. Development and evaluation of soil moisture deficit index and evapotranspiration deficit index for agricultural drought monitoring. Agricultural and Forest Meteorology 133: 69 88. Nearing, M. A., L. Deer Ascough, and J. M. Laflen , 1990. Sensitivity analysis of the WEPP hillslo pe profile erosion model. Transactions of ASAE 33(3): 839 849. Neitsch, S. L., J. G. Arnold, J.R. Kiniry, and J. R. Williams, 2002. Soil and Water Assessment Tool user's manual, Version 2000. Temple, Tex.: USDA ARS Grassland, Soil, and Water Research Laboratory. Available at: www.brc.tamus.edu/swat/ . (accessed in n July 2012) Nossent, J., P. Elsen, and W. Bauwens, 2011. environmental model. Environmental Modelling and Software 26: 1515 1525. Oakley, J.E. and 2004. Probabilistic sensitivity analysis of complex models: A Bayesian approach. Journal of the Royal Statistical Society Series B 66: 751 769. doi: 10.1111/j.1467 9868.2004.05304.x Osidele, O.O., W. Zheng, and M.B. Beck, 2003. Coping with uncertainty : A case study in sediment transport and nutrient load analysis. Journal of Water Resources Planning and Management 4 : 345 355. Padilla, F.L.M., M.P. Gonzalez Dug o, P. Gavilan, and J. Dominguez, 2011. Integration of vegetation indices into a water balance model to estimate evapotranspiration of wheat and corn. Hydrology and Earth System Science 15: 1213 1225. Palmer, W. C . , 1968. Keeping track of moisture conditions, nationwide: the new crop moisture index. Weatherwise 21: 156 161. Palmer, W. C. , 1965. Meteorological drought. Research Paper No. 45. Washington, D. C.: U. S. Department of Commerce. PAGE 199 199 Pappenberger, F., K. J. Beven, M. Ratto, and P. Matgen, 2008. Multi method global sensitivity analysis of flood inun dation models. Advances in Water Resources 31: 1 14. Pappenberger, F., I. Iorgulescu, and K. J. Beven, 2006. Sensitivity analysis based on regional splits and regre ssion trees (SARS RT). Environmental Modelling and Software 21(7): 976 990. Passioura, J. B. , 1983. Roots and drought resistance. Agricultural Water Management 7: 265 280. Patil, A ., Z.Q. Deng, and R. F. Malone, 2011. Input data resolution induced uncertaint y in watershed m odeling. Hydrological Processes 25: 2302 2312. Prats, A.G. and S.G. PicÃ³, 2010. Performance evaluation and uncertainty measurement in irrigation scheduling soil water balance approach. Journal of Irrigation and Drainage Engineering ASCE doi : 10.106/(ASCE)IR.1943 4774.0000245. R, version 1.9. http://cran.r project.org/web/packages/sensitivity/sensitivity.pdf (a cc essed Se ptember 2014). Pujol, G. , 2009. Simplex based screening designs for estimating metamodels. Reliabilit y Engineering and System Safety 94: 1156 1160. Raes, D., P. Sted uto, T.C. Hsiao, and E. Fereres, 2009. AquaCrop The FAO crop model to simulate yield response to water. II. Main algorithms and software description. Agronomy Journal 101: 438 447. Ratliff, L. F. , J.T. Ritchie, and D.K. Cassel, 1983. Field measured limits of soil water availability as related to laboratory measured prope rties. Soil Science Society of America Journal 47: 770 775. Ratto, M., S. Tarantola, and A. Saltelli , 2001. Sensitivity analysis in model calibration: GSA GLUE a pproach. Computer Physics Communications 136: 212 224. Rawls, W.J. a nd D.L. Brakensiek, 1989. Estimation of soil water retention and hydraulic properties, in unsaturated flow in hydrologic modeling: Theory and Practice, edited by H. Morel Seytoux, 275 300, Kluwer Acad, Dordrecht. Refsgaard, J.C. and B. Storm, 2012. Chapter 23 :MIKE SHE. In Computer Models of Watershed Hydrology. Ed V.P. Singh. Water Resources Publications, Highland Ranch, CO. Refsgaard, J.C., J.P. van der Sluijs, A.L. HÃ¸jberg, P.A. Vanrolle ghem, 2007. Uncertainty in environmental modeling process A framework and guidance. Environmental Modelling and Software : 1543 1556. Restrepo, P. and J. Schaake , 2009. Hydrologic forecasting in 21 st century: Challenges and directions of research, Europe an Geophysical Union Assembly, vol 11. EGU2009 5574 1 PAGE 200 200 Reusser, D.E., W. Buytaert, and E. Zehe, 2011. Temporal dynamics of model parameter sensitivity for computationally expensive models with the Fourier amplitude sensitivity test. Water Resource Research 47. Ritchie, J. T., A. Gerakis, and A. Suleiman , 1999. Simple model to estimate field measured soil water limit, Transactions of ASAE 42: 1609 1614. Carpena, 2013. Performance evaluation of hydrologic model: Statistical significanc e for reducing subjectivity in goodness of fit assessment. Journal of Hydrology 480: 33 45. Robertson, M.J., S. Fuka i, M.M. Ludlow, and G.L. Hammer, 1993a. Water extraction by grain sorghum in a sub humid environment. I. Analysis of the water extraction pa ttern. Field Crops Research 33: 81 97. Robertson, M.J., S. Fuka i, M.M. Ludlow, and G.L. Hammer, 1993b. Water extraction by grain sorghum in a sub humid environment. II. Extraction in relation to root growth. Field Crops Research 33: 99 112. Romero, C.C. an d M.D. Dukes, 2009. Turfgrass and ornamental plant evapotranspiration and crop coefficient literature review, Agricultural and Biological Engineering Department, University of Florida, Gainesville, Florida. http://abe.ufl.edu/mdukes/pdf/turfgrass/Turfgrass Final Review.pdf (Accessed on 15 December 2012). Ross, M., J. Geurink, A. Said, A. Aly, and P. Tara, 2005. Evapotranspiration conceptualization in HSPF MODFLOW integrated models. Journal of the AWRA 41(5): 1013 1025. Ruano, M.V, J. Ribes, A. Seco, and J. Ferrer , 2012. An improved sampling strategy based on trajectory design for application of the Morris method to systems with many input factors. Environmental Modelling and Software 37: 103 109. Salemi, A., M.A. Soom, T.S. Lee, S.F. Mousav i, A. Ganji, and M. KamilYusoff, 2011. Application of AquaCrop model in deficit irrigation management of winter wheat in Arid region. African Journal of Agricultural Research 610: 2204 2215. Saltelli, A. and P. Annoni, 2010. How to avoid a perfunctory sensitivity analysis. Environmental Modelling and Software 25: 1508 1517. doi: 10.1016/j.envsoft.2010.04.012 Saltelli, A., P. Annoni, I. Azzini, F. Campolo ngo, M. Ratto, and S. Tarantola, 2010. Variance based sensitivity analysis of model output. Design and estimator for total sensitivity index. Computer Physics Communications 181: 259 270. Saltelli, A., F. Campolongo, and J. Cariboni, 2009. Screening important inputs in models with strong intera ction properties. Reliability Engineering and System Safety 94: 1149 1155. doi:10.1016/j.ress.2008.10.007 PAGE 201 201 Saltelli, A., M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Galtelli, M. Saisana, and S. Tarantola , 2008 . Global Sensitivity Analysis The Primer , Chichester, U.K., John Wiley and Sons, ISBN 978047005997 . Saltelli, A., M. Ratto, S. Tarantola, and F. Campolongo , 2005 . Sensitivity analysis for chemical models . Chemical Rev iews 105(7): 2811 2827 . Saltelli, A., S. Taranto la, F. Campolongo, and M. Ratto , 2004 . Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models . Chichester, U.K.: John Wiley and Sons, ISBN 9780470870938 . Saltelli, A. , 2002. Making best use of model valuations to compute sensitivity indices. Computer Physics Communications 145: 280 297. Saltelli, A., K. Chan, and E.M. Scott , 2000 . Sensitivity Analysis, Chichester U.K., John Wiley and Sons, ISBN: 0471998923 . Sal telli, A. , S. Tarantoal, and K. P. S. Chan , 1999. A quantitative model independent method for global sensitivity analysis of model outpu t. American Society for Quality 41(1): 39 56. Saltelli, A. and 1995. About the use of rank transformation in sensitivity analysis of model output. Reliability Engineering and System Safety 50: 225 239. Saltelli, A. and J. Marivoet , 1990. Non parametric statistics in sensitivity analysis for model output: A comparison of selected techniques . Reliability Engineering and System Safety 28: 229 253. Santer, T.J., B.J. William, and W.I. Notz, 2003. The design and analysis of computer experiments. Berlin: Springler. Santiago, J., M. Cla eys Bruno, and M. Sergent, 2012a. Const ruction of spacefilling designs using WSP algorithm for high dimensional spaces. Chemometrics and Intelligent Laboratory Systems 113: 26 31. doi: 10.1016/j.chemolab.2011.06.003 Santiago, J., B. Corre, M. Claeys Bru no, and M. Sergent, 2012b. Improved sensit ivity through Morris extension. Chemometrics and Intelligent Laboratory Systems 113:52 57. doi: 10.1016/j.chemolab.201110.006 Satti, S.R., J.M. Jaco bs, and S. Irmak, 2004. Agricultural water management in a humid region: sensitivity to climate, soil and cr op parameters. Agricultural Water Management 70: 51 65. Saxton, K.E. and J.W. Rawls, 2006. Soil water characteristics estimates by texture and organic matter for hydrologic solutions. Soil Science Society of America Journal 62(4): 847 855. Saxton, K.E., W.J. Rawls , J.S. Romberger, and Papendick, 1986. Estimating generalized soil water characteristics from texture. Transactions of ASABE 50: 1031 1035. PAGE 202 202 Shen, Z., Q. Hong, H. Yu, and R. Liu, 2008. Parameter uncertainty analysis of the non point source polluti on in the Danning River watershed of the Three Gorges Reservoir Region, China. Science of Total Environment 405: 195 205. Shin, M J., J.H. A. Guillaume, B.F. W. Croke, and A. J. Jakeman, 2013. Addressing ten questions about conceptual rainfall runoff models w ith global sensitivity analysis in R. Journal of Hydrology 503: 135 152. z Carpena, C. Dharmasri, A. Saxton, M. Arabi, M.L. Wolfe, J. Frankenberger, C. Graff , and T.M. Sohrabi , 2006. Uncertai nty in TMDL models. Transactions of ASABE 49(4): 1033 1049. Sie ber, A. and S. Uhlenbrook, 2005. Sensitivity analyses of a distributed catchment model to verify the model structure. Journal of Hydrology 310(1 4): 216 235. SimLab. 2008. Version 2.2 Reference Manual European Commission, Joint Research Center, Ispra. http://sensitivity analysis.jrc.cec.eu.int/ (Assessed on 10 December 2011). Singh, V.P., 2012. Computer Models of Watershed Hydrology. Water Resources Publication. Highland Ranch, C O , ISBN: 1887201742. Singh, V. P. and D. K. Frevert , 2006. Watershed Models . Taylor and Francis Group, Boca Raton, F L , ISBN 0 8493 3609 0. Singh, V. P. and D.K. Frevert , 2002a. Mathematical Models of Large Watersh ed Hydrology . Water Resources Publications , Highlands Ranch, C O . Singh, V. P. and D.K. Frevert , 2002b. Mathematical Models of Small Watershed Hydrology and Applications . Water Resources Publications, Highlands Ranch, C O . Smith, R.M.S. and H.S. Wheater, 2004. Multiple objective evaluation of a simple phosphorus transfer model. Hydrological Processes 18: 1703 1720. and S. Kucherenko, 2009. Derivative based global sensitivity measures and their link with global sensitivity indices. Mathematics and Computers in Simulation 79: 3009 3017. , 1993. Sensitivity analysis for non l inear mathematical models. Mathe matical Modelling and Comput er Experiments 1: 407 414. , 1976. Uniformly distributed sequences with an additional uniform property. USSR Computational Mathematics and Mathematical Physics 16: 236 242. Sohier, H., J . L. Farges, and H. Piet Lahan ier, 2014. Improvement of the representativity of the Morris method for air launch to orbit separation. Reprints of the 19 th World Cong ress, Cape Town, South Africa. PAGE 203 203 Soil Survey Staff. 2004. National soil characterization database. USDA/NRCS National Soil Survey Center, Lincoln, NE. Soil and Water Engineering Technology, Inc. (SWET), 2014. Watershed Assessment Model Model Documentation. SWET, 2011. Watershed Assessment Model Documentation and Validation: Model Calibration and Validation Report. Submitted to SFWMD. SWET, 2009 a . WAM enhancement and application in the Lake Okeechobee Watershed. Preparared in association with HDR. SWET, 2009b. Upper Peace River Three Lakes Modeling Report using WAM and WASP. Complete for Tetra Tech Inc. SWET, 2008. Watershed Assessment Model Model Documentation. SWET, 1997. Development of a GRID GIS based simulation model for the lower St. Johns River Basin Deep Creek sub basin hydrologic/water quality assessment Final Report. Submitted to St. Johns River Water Management District. Sorooshian, S. and H.V. Gupta , 1995. Model calibration. Chapter 2 in Computer models of Watershed Hydrology, edited by Singh, Water Resources Publication . South Florida Water Management District (SFWMD), 2010. Nutrient budget analysis for the Lake Okeechobee watershed. Prepared in association with HDR, UF IFAS, JGH Engineering, and SWET Inc. SFWMD, 2009. Development of design criteria for stormwater treatment areas (STAs) in the Northern Lake Okeechobee Watershed, Final Report. Spear , R.C. and G .M. Hornberger, 1980. Eutrophication in Peel Inlet, II, Identification of critical uncertainties via Generalised Sensitivity Analysis. Water Research 14: 43 49 . Sr Carpena, and R.M. Maxwell, 2014. Insights on geologic an d vegetative control over hydrologic behavior of a large complex basin Global Sensitivity Analysis of an integrated parallel hydrologic model. Journal of Hydrology (accepted) . Srivastava, V., 2013. Geologic, vegetative and climatic controls on coupled hy drologic proesses in a complex river basin: Lessons learned from a fully integrated hydrologic model. Ph.D. dissertation. University of Florida, Gainesville, FL. Srivastava, P., J.M. Harmel, P.D. Robbilard , and R.L. Day , 2002. Watershed optimization of bes t management practices using A nn AGNPS and a generic algorithm . Water Resource Research 38(3): 1 14 . PAGE 204 204 Steduto, P., T.C. Hsiao, D. Raes, and E. Fereres, 2009. AquaCrop The FAO crop model to simulate yield response to water. I. Concepts and underlying principles. Agronomy Journal 101: 426 437. Stow C.A., K.H. Reckhow, S.S. Qian, E.C. Lamon III, G.B. Arhonditsis, M.E. Borsuk, and D. Seo , 2007. Approaches to evaluate water quality model parameter uncertainty for adaptive TMDL implementation. Journal of the AWRA 43(6): 1499 1507. Sudheer, K. P., I. Chaubey , V. Garg, and K. W. Migliaccio, 2007. Impact of time scale of the calibration objective function on the per forma nce of watershed models. Hydrological Processes 21: 3409 3419. Suleiman, A. and J. T. Ritchie, 2004. Modifications to the DSSAT vertical drainage model for more accurate soil water dynamics estimation. Soil Science 169 (11): 745 757. Tang, Y., P. Reed, K. van Werkhoven, and T. Wagener , 2007. Advancing the identification and evaluation of distributed rainfall runoff models using global sensitivity analysi s. Water Resources Research 43: W06415, doi:10.1029/2006WR005813 Tarantola, S. and A. Saltelli. 2003. SAMO 2001: Methodological advances and innovative applications of sensitivity analysis. Reliability Engineering and Systems Safety 79: 121 122. Thornthwaite, C. W., 1948. An approach toward a rational classification of climate. Geographical Review 38: 55 9 4. Thornton, P.E. and S.W. Running, 1999. An improved algorithm for estimating incident daily solar radiation from measurements of temperature, humidity, and precipitation. Agricultural and Forest Meteorology 93: 211 228. Thornton, P.E ., S.W. Running, and M.A. White, 1997. Generating surfaces of daily meteorology variables over large regions of complex terrain. Journal of Hydrology 190: 214 251. Tong, C.H. and F.R. Graziani, 2007. A global sensitivity analysis methodology for multi physics applications. Law rence Livermore National Laboratory, CA, USA, report number UCRL TR 227800. United States Department of Agriculture, 1986. Urban Hydrology for Small Watersheds. Technical Release 55 (TR 55). United States Environmental Protection Agency (US EPA), 2013. Eco system research: Watershed and wate quality modeling technical center, Athens, GA. http://www.epa.gov/athens/wwqtsc/index.html (accessed October 2013). US EPA. 2010. Water Quality Standards for th e State of Florida's Lakes and Flowing Waters. Fact Sheet. http://water.epa.gov/lawsregs/rulesregs/florida_factsheet.cfm#summary (accessed September 2012). PAGE 205 205 US EPA, 2009 a . Guidance on the development, evaluation, and application of environmental models, Technical Report, Office of Science Advisor, Council for Regulatory Environmental Modeling. EPA/100/K 09/003. http://www.epa.gov/crem/library/cred_guidance_0309.pdf. (accessed November 2013) . US EPA, 2009 b . Proposed Total Maximum Daily Loads (TMDLs) for the Lake Alfred, Crystal Lake, Lake Arianna North nutrients. September 2009. United States Geologi cal Survey (USGS), 2014. National Water Data for the Nation. http://waterdata.usgs.gov/nwis (accessed August 2014). v an Griensven, A., T. Meixner, S. Grunwald, T. Bishop, M. Diluzio, and R. Srinivasan, 2006. A global sensitivity analysis tool for the parameters of multi variable catchment models. Journal of Hydrology 324(1 4): 10 23. van Werkhoven, K., T. Wagener, P. Reed, and Y. Tang, 2009. Sensitivity guided reduction of parametric dimensionality for multi ob jective calibration of watershed models. Advances in Water Resources 32(8): 1154 1169. doi: 10.1016/j.advwatres.2009.03.002 van Werkhoven, K., T. Wagener, P. Reed, and Y. Tang, 2008. Rainfall characteristics define the value of streamflow observations for distributed watershe d model identification. Geophyscial Research Letters 35: L 11403. doi:10.1029/2008GL034162 Vezzaro, L. and P .S. Mikkelsen, 2012. Application of global sensitivity analysis and uncertainty quantification in dynamic modelling of micro pollutant in stormwater runoff. Environmental Modelling and Software 27 28: 40 51. Vogel, R.M., 2006. Regional calibration of water shed models. In Watershed Models Ed. Singh and Frevert, ISBN 0 8493 3609 0. Vrugt, J.A. C.J.F. er Braak, H.V. Gupta, and B.A. Robinson, 2009. Equifinality of formal (DREAM) and informal (GLUE) Bayesian approaches in hydrologic modeling? Stochastic Environm ental Resources Risk Assessment 23:101 1026. doi:10.1007/s00477 008 0274 y. Vrugt, J.A., C.J.F. ter Braak, M.P. Clark, J.M. Hyman, and B. A. Robinson, 2008. Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Morkov Chain Mo nte Carlo simulation. Water Resources Res earch 44: W00B09. doi: 10.1029/2007WR006720. Wagener, T., K. van Werkhoven, P. Reed, and Y. Tang, 2009. Multiobjective sensitivity analysis to understand the information content in streamflow observations for distributed watershed modeling. Water Resources Research 45: W02501. doi:10.1029/2008WR007347. Wagene r, T. and J. Kollat , 2007. Numerical and visual evaluation of hydrological and environmental models using Monte Carlo analysis toolbox. Environmental Modeling and S oftware 22: 1021 1033. PAGE 206 20 6 Wagener, T., H.W. Wheater , and H.V. Gupta , 2004 . Rainfall runoff mode ling in gages and ungagged catchements, Imperial College Press. Wainwright, H.M., S. Finsterle, Y. Jung, Q. Zhou, and J.T. Birkholzer, 2014. Making sense of global sensitivity anal yses. Computers and Geosciences 65: 84 94. Wallace, S.D., and R.L. Knight, 2006. Small scale constructed wetland treatment systems feasibility, design criteria, oeration , and maitainance requirement, Final Report. Wallach, D., D. Mako wski, J.W. Jones, and F. Brun, 2014. W orking with Dynamic Crop Models, Second Edition. Academic Press, Elsevier. Wallach , D., D. Makowski, and J.W. Jones, 2006. W orking with Dynamic Crop Models, First Edition . Amsterdam, The Netherlands: Elsevier . Wang, X. , S.R. Potter, J.R. Will iams, J.D. Atwood, and T. Pitts, 2006. Sensitivity analysis of APEX for natio nal assessment, Transactions of ASABE 49(3): 679 688. Warmink, J.J., J.A.E.B. Janssen, M.J. Booij, and M.S. Krol, 2010. Identification and classification o f uncertainties in application of environmental models. Environmental Modelling and Software 25: 1518 1527. doi:10.1016/j.envsoft.2010.04.011 . Wiens, D.P. , 1991. Designs for approximately linear regression: Two optimality properties of uniform designs. Statistics and Probability Letters 12: 217 221. Wilcoxon, F. , 1945. Individual comparisons by ranking methods. Biometrics 1: 80 83. White, K. L. and I. Chaubey , 2005. Sensitivity analysis, calibration, and validations for a multisite a nd multivariable SWAT model. Journal of the AWRA 41(5): 1077 1089. Whiting, W.B., T.M. Tong, and M.E. Ree d, 1993. Effect of uncertainties in thermodynamic data model parameters on calculated process performance. Industrial and Engineering Chemistry Research 32: 1367 1371. Woli, P., J.W. Jones, and K.T. Ingram, 2013. Assessing the Agricultural Reference Index for Drought (ARID) using uncertainty and sensitivity analyses. Agronomy Jou rnal 105(1): 150 160. Woli, P., J.W. Jones, K.T. Ingram, and C.W. Fraisse, 2012. An Agricultural Reference Index for Drought. Agronomy Journal 104 (2): 287 300. WÃ¶sten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas, 1999. Development and use of a database of hydraulic properties of European soils. Geoderma 90(3 4): 169 185. sensitivity analysis and its application to tree g rowth modeling. Reliability Engineering and System Safety 107: 35 43. PAGE 207 207 Wu, Y. and S. Liu; 2012. Auomatic calibration, sensitivity and uncertainty analyses of a complex model using the R package Flexible Modeling Environment (FME): SWAT as an ex ample Enviro n mental. Modelling and Software 31: 99 101 . Wu, S., J. Li, and G. H. Huang, 2007. Characterization and evaluation of elevation data uncertainty in water resource modelin g with GIS. Water Resources Research 22: 959 972. Wu, S., J. Li, and G. H. Huang, 2005. An evaluation of grid size uncertainty in empirical soil loss modeling with di gital elevation models. Environmental Modelling and Softw are 10: 33 42. Wu, J., R. Zhang, and S. Gui, 1999. Modeling soil water movement with water uptake by roots. Plant and Soil 215: 7 17. Wu, L. , 1985. Matching irrigation to turfgrass root depth, California turfgrass culture, University of California, 35: 1 4. Xu, C. and G.Z. Gertner, 2008. Uncertainty and sensitivity analysis for models with correlated parameters. Reliability Engineering and System Safety 93: 1563 1573. Yang, J ., Y. Liu, W. Yang, and Y. Chen, 2012. Multi objective sensitivity analysis of a fully distributed model Wet Spa. Water Reseourc e Management 26: 109 128. Yu, P. S., T. C. Yang, and S. J. Chen, 2001. Comparison of uncertainty analysis methods for a distributed rainfall runoff model. Journal of Hydrology 244: 43 59. Yuan, Y., Y.P. Khare, X. Wang, X., P. Parajuli, I. Kisseka, and S. Finsterle, 2014. Hydrologic and water quality mode ls: Sensitivity analysis. Transactions of ASABE (in revision). Zajac, Z., 2010. Global sensitivity and uncertainty analysis of spatially distributed watershed model. Ph.D. Dissertation. Dept. of Agri. And B io. Eng., University of Florida, Gainesville, FL. Zhan, C S., X. M. Song, J. Xia, and C. Tong, 2013. An efficient integrated approach for global sensitivity analysis of hydrological model parameters. Environmental Modelling and Software 41: 39 52. doi.org/ 10.1016/j.envsoft.2012.10.009 Zhang, J., and B. Sharfstein, 2013. Chapter 8: Lake Okeechobee Watershed protection program. Zhang, J., J.G. Hiscock, A.B. Bottcher, B.M. Jacobson, P.J. Bohlen, 2006. Modeling Phosphorus load reduction of agricultural water management practice on a beef cattle ranch. ASABE Annual International Meet Paper # 062010, Portland, OR. Zhang, J. and S.I. Gornak, 1999. Evaluation of field scale water quality models for the lake Okeec hobee regulatory program. Transactions of ASAE 15(5): 441 447. PAGE 208 208 Zhang, K., R. Srinivasan, and D. Bosch, 2009 . Calibration and uncertainty analysis of the SWAT model using Genetic Algorithms and Bayesian Model Averaging. Journal of Hydrol ogy 374: 307 317. Zheng, Y. and A. A. Keller. 2006. Understanding parameter sensitivity and its management implications in watershed scale water quality modeling. Water Resources Research 42(1): 1 14. PAGE 209 209 BIOGRAPHICAL SKETCH Yogesh Khare was born an d brough t up in Mumbai, India. He received his B achelor of Civil Engineering from the University of Mumbai in 2006. After working in Larsen and Toubro ECC Division as a trainee & senior engineer and Sardar Patel College of Engineering as a lecturer he came to the University of Florida (UF) to pursue his Master of Science in Civil and Coastal Engineering. He started working in the Department of Agricultural and Biological Engineering (ABE) at UF as a part time research assistant in 2009. After receiv ing mast in Coastal Engineering in early 2010 he became a full time doctoral student in the Hydrologic Modeling group at ABE . |