Water Balance Analysis Using WARMF in Support of an EFDC Hydrodynamic Model by Thomas Weems A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master?s of Science in Civil Engineering Auburn, Alabama August 3, 2013 Keywords: water balance, WAMRF, hydrology, CE-QUAL-W2, EFDC, watershed Copyright 2013 by Thomas Weems Approved by Xing Fang, Chair, Associate Professor of Civil Engineering Jose G. Vasconcelos, Assistant Professor of Civil Engineering Luke J. Marzen, Professor of Geology and Geography ii Abstract An Environmental Fluid Dynamics Code (EFDC) model was previously developed for Bankhead Reservoir, located in northwest Alabama. The EFDC model inaccurately predicted water surface elevations in the reservoir, due to sources of mass that were not able to be quantified. The reservoir has numerous tributaries contributing flow to the domain, as well as a known leakage issue through the spill gates at the downstream Bankhead Dam. Because EFDC does not include a mass balance tool, the objective of the study is to use other modeling tools (hydrologic modeling and water balance analysis) to account for the unknown mass inputs to the reservoir. The watershed and tributary inflows were modeled using Watershed Analysis Risk Management Framework (WARMF). A traditional hydrologic modeling approach to optimize the modeling of peak flows was used initially, but found to be inadequate for the purpose of predicting flows for the specific low-flow time periods in the EFDC model. Therefore, separate 2010 and 2011 low-flow WARMF models were developed to achieve the best possible re-creation of observed data for the time periods of interest. Finally, to bring the reservoir into balance and to quantify downstream leakage, the water balance utility contained within CE-QUAL-W2 was coupled to the EFDC model to determine a flow file that would bring the reservoir into balance. After the balancing procedure was completed, the EFDC model simulated the water surface elevations to a high accuracy, with absolute errors less than 15 cm. iii Acknowledgments I would like to thank my advisor, Dr. Xing Fang, for the opportunity to conduct research under his guidance. This project would not have been possible without his instruction, patience, and thorough knowledge of the applicable fields. His ability to instruct and motivate made graduate study as enjoyable as it was instructional. I would also like to thank my other committee members, Dr. Jos? Vasconcelos and Dr. Luke Marzen, for their gifted instruction which prepared me to undertake this project. Also, this project would have been extremely difficult without the willing support of the other students at Auburn, Janesh Devkota, Gang Chen, Manoj KC, and Liping Jiang. I would also like to thank Dr. Bill Garrett, Jon Ponstein, and my other co-workers for their support and also for providing valuable data. To my family, friends, and fianc?, I appreciate your vast patience and support while I was working on this project, even though I might not have always returned in kind. I am blessed to be surrounded by such a good support system at home, at school, and at work. iv Table of Contents Abstract ................................................................................................................................ ii Acknowledgments ............................................................................................................... iii List of Tables .................................................................................................................... viii List of Figures ..................................................................................................................... ix List of Abbreviations ......................................................................................................... xiii Chapter 1 Introduction and Project Purpose ...........................................................................1 1.1 Background ................................................................................................................1 1.2 Scope and Objectives ..................................................................................................1 1.3 Study Area ..................................................................................................................4 1.4 Thesis Organization ....................................................................................................7 Chapter 2 Development of WARMF Model ..........................................................................8 2.1 Introduction ................................................................................................................8 2.2 Model Description ......................................................................................................9 2.3 Literature Review ..................................................................................................... 13 2.4 Methods .................................................................................................................... 19 2.4.1 Watershed Delineation Using Basins .................................................................. 19 v 2.4.2 WARMF Model Input ........................................................................................ 24 2.4.2.1 Observed Streamflow Data .............................................................................. 25 2.4.2.2 Atmospheric Data ........................................................................................... 25 2.4.2.3 Soil Data ......................................................................................................... 28 2.4.3 WARMF Simulation .......................................................................................... 29 2.4.4 Calibration Criteria ............................................................................................ 30 2.4.4.1 Literature Review of Hydrologic Model Calibration ........................................ 31 2.4.4.2 Coefficient of Determination (R2) ................................................................... 31 2.4.4.3 Nash-Sutcliffe Efficiency (NSE) ..................................................................... 32 2.4.4.4 RMSE and Standard Deviation Ratio (RSR) .................................................... 33 2.4.4.5 Percent Bias (PBIAS) ...................................................................................... 34 2.4.4.6 Evaluation Statistics Used ............................................................................... 34 2.5 Results and Discussion ............................................................................................. 36 2.5.1 Statistical Summary ........................................................................................... 37 2.5.2 Blackwater Creek ............................................................................................... 38 2.5.3 Locust Fork ........................................................................................................ 43 2.5.4 Lost Creek ......................................................................................................... 47 2.5.5 Mulberry Fork .................................................................................................... 51 2.5.6 Discussion of Calibration and Validation Results ............................................... 55 2.6 Summary .................................................................................................................. 59 vi Chapter 3 Water Balance Analysis....................................................................................... 61 3.1 Introduction/Background .......................................................................................... 61 3.2 Water Balance Development Using WARMF 2010 and 2011 Low-flow Models ...... 62 3.2.1 Development of the WARMF Low-Flow Hydrologic Model .............................. 65 3.2.2 WARMF Low-flow Model Setup ........................................................................... 68 3.2.2.1 Recalibration of WARMF Model for Low-flow Conditions ............................ 68 3.2.2.2 Analysis of 2010 and 2011 Low-flow Calibrations .......................................... 68 3.3 Background of CE-QUAL and EFDC Hydrodynamic Models .................................. 80 3.3.1 CE-QUAL-W2 ................................................................................................... 80 3.3.2 CE-QUAL-W2 Water Balance Utility ................................................................ 81 3.3.3 EFDC (Environmental Fluid Dynamics Code) Hydrodynamic Model ................ 84 3.4 Development of EFDC Model to Incorporate Water Balance .................................... 85 3.4.1 EFDC Boundary Condition Setup ...................................................................... 87 3.5 EFDC Model Methodology ....................................................................................... 98 3.6 EFDC Model Results and Analysis ......................................................................... 101 3.6.1 2010 BALANCED EFDC Model ...................................................................... 101 3.6.2 2011 BALANCED EFDC Model ..................................................................... 105 3.6.3 Summary of EFDC Modeled Elevations ........................................................... 106 Chapter 4 Summary and Future Study ............................................................................... 116 4.1 Summary of Study .................................................................................................. 116 vii 4.2 Recommendations for Future Study ........................................................................ 118 References ......................................................................................................................... 119 viii List of Tables Table 2.1 Land Cover data from sample catchments ........................................................... 23 Table 2.2 USGS gage stations located in Bankhead watershed ............................................ 25 Table 2.3 Soil data downloaded from SSURGO for the three major counties ...................... 29 Table 2.4 Statistical benchmarks for evaluating hydrologic models (Moriasi et al. 2007) .... 36 Table 2.5 Statistical summary of calibration and validation for streamflow ......................... 37 Table 2.6 Calibration coefficients for the Blackwater Creek subwatershed ......................... 40 Table 2.7 Calibration coefficients for the Locust Fork subwatershed .................................. 47 Table 2.8 Calibration Coefficients for the Lost Creek Subwatershed ................................... 51 Table 2.9 Calibration coefficients for the Mulberry Fork subwatershed .............................. 55 Table 2.10 Calibrated Parameters Compared to Literature Values....................................... 57 Table 3.1 Calibration Parameters for WARMF Low-flow Calibration ................................ 68 Table 3.2 Statistical results for the 2010 and 2011 Low-flow models .................................. 69 Table 3.3 Statistical summary of elevation differences between the 2010 NO BALANCE and 2010 BALANCED models .............................................................................. 105 ix List of Figures Fig. 1.1 Overview of Bankhead Reservoir in relation to neighboring counties .......................5 Fig. 1.2 Location of available meteorological data in towns near Bankhead Reservoir ..........5 Fig. 1.3 Tributaries contributing to the Black Warrior River and Bankhead Dam .................6 Fig. 1.4 Locations of USGS gage stations in the Bankhead watershed ...................................6 Fig. 2.1 Hydrologic components within WARMF (Chen et al. 2001) .................................. 10 Fig. 2.2 Representation of hydrologic cycle (Dardashti 2010) .............................................. 14 Fig. 2.3 HUC units contributing to Bankhead Watershed .................................................... 20 Fig. 2.4 Average slopes of delineated catchments ............................................................... 21 Fig. 2.5 Layout of delineated catchments with river segments and reservoir in WARMF .... 22 Fig. 2.6 Location of sample catchments for land cover data ................................................ 24 Fig. 2.7 Location of meteorological stations in watershed region ........................................ 27 Fig. 2.8 Blackwater Creek subwatershed ............................................................................ 38 Fig. 2.9 Flow calibration of Blackwater Creek subwatershed .............................................. 41 Fig. 2.10 Flow validation of Blackwater Creek ................................................................... 42 Fig. 2.11 Locust Fork subwatershed ................................................................................... 44 Fig. 2.12 Flow calibration of Locust Fork subwatershed ..................................................... 45 Fig. 2.13 Flow validation of Locust Fork subwatershed ...................................................... 46 Fig. 2.14 Lost Creek subwatershed ..................................................................................... 48 Fig. 2.15 Calibration of Lost Creek subwatershed ............................................................... 49 x Fig. 2.16 Validation of Lost Creek subwatershed ................................................................ 50 Fig. 2.17 Mulberry Fork subwatershed ............................................................................... 52 Fig. 2.18 Calibration of Mulberry Fork subwatershed ......................................................... 53 Fig. 2.19 Validation of Mulberry Fork subwatershed .......................................................... 54 Fig. 3.1 2010 High-flow calibration of Blackwater Creek over low-flow period ................. 63 Fig. 3.2 2010 High-flow calibration of Locust Fork over low-flow period .......................... 63 Fig. 3.3 2010 High-flow calibration of Lost Creek over low-flow period ............................ 64 Fig. 3.4 2010 High-flow calibration of Mulberry Fork over low-flow period ...................... 64 Fig. 3.5 Sum of the four USGS tributaries compared to the total USGS observed flows in 2010 ........................................................................................................................ 65 Fig. 3.6 Calibration of Blackwater Creek for 2010 low-flow model .................................... 72 Fig. 3.7 Calibration of Blackwater Creek for 2011 low-flow model .................................... 73 Fig. 3.8 Calibration of Locust Fork for 2010 low-flow model ............................................. 74 Fig. 3.9 Calibration of Locust Fork for 2011 low-flow model ............................................. 75 Fig. 3.10 Calibration of Lost Creek for 2010 low-flow model ............................................. 76 Fig. 3.11 Calibration of Lost Creek for 2011 low-flow model ............................................. 77 Fig. 3.12 Calibration of Mulberry Fork for 2010 low-flow model ....................................... 78 Fig. 3.13 Calibration of Mulberry Fork for 2011 low-flow model ....................................... 79 Fig. 3.14 Comparison of EFDC and CE-QUAL bathymetries for the whole reservoir ......... 83 Fig. 3.15 View of CE-QUAL model near Bankhead Dam ................................................... 83 Fig. 3.16 View of CE-QUAL model near Plant Gorgas ....................................................... 84 Fig. 3.17 Comparison of 2010 WARMF modeled flows of Locust Fork at the USGS gage and at the intersection with Bankhead Reservoir ...................................................... 86 xi Fig. 3.18 Locations of boundary condition tributaries in the EFDC model domain. Major locations of interest are labeled. Other points are minor tributaries or the pour point of a subcatchment located on the reservoir. .................................................................. 87 Fig. 3.19 Upstream boundary of EFDC model .................................................................... 89 Fig. 3.20 Downstream boundary conditions of EFDC model .............................................. 89 Fig. 3.21 2010 upstream flow boundary at Smith Dam ....................................................... 90 Fig. 3.22 2011 upstream flow boundary at Smith Dam ....................................................... 90 Fig. 3.23 2010 downstream flow boundary at Bankhead Dam............................................. 91 Fig. 3.24 2011 downstream flow boundary at Bankhead Dam............................................. 91 Fig. 3.25 2010 air temperature from Birmingham Airport ................................................... 92 Fig. 3.26 2011 air temperature from Birmingham Airport ................................................... 92 Fig. 3.27 2010 relative humidity from Birmingham Airport ................................................ 93 Fig. 3.28 2011 relative humidity from Birmingham Airport ................................................ 93 Fig. 3.29 2010 rainfall data from Birmingham Airport ........................................................ 94 Fig. 3.30 2011 rainfall data from Birmingham Airport ........................................................ 94 Fig. 3.31 2010 solar radiation data from Birmingham Airport ............................................. 95 Fig. 3.32 2011 solar radiation data from Birmingham Airport ............................................. 95 Fig. 3.33 2010 wind speed data from Birmingham Airport ................................................. 96 Fig. 3.34 2011 wind speed data from Birmingham Airport ................................................. 96 Fig. 3.35 Withdrawal/return boundary conditions simulating Plant Gorgas ......................... 97 Fig. 3.36 Intake and discharge paths for Plant Gorgas ......................................................... 97 Fig. 3.37 2010 Mass balance flow boundary added near Bankhead Dam .......................... 100 Fig. 3.38 2011 Mass balance flow boundary near Bankhead Dam ..................................... 100 Fig. 3.39 Location of Cordova USGS gage ....................................................................... 101 xii Fig. 3.40 2010 comparison of modeled and observed water surface elevations at Bankhead Dam ....................................................................................................................... 108 Fig. 3.41 2010 comparison of final balanced water surface elevation to the observed elevation at Bankhead Dam.................................................................................... 109 Fig. 3.42 2010 comparison of modeled and observed water surface elevations at Cordova USGS station ......................................................................................................... 110 Fig. 3.43 2010 comparison of final balanced water surface elevation to the observed elevations at Cordova USGS station ....................................................................... 111 Fig. 3.44 2011 comparison of modeled and observed water surface elevations at Bankhead Dam ....................................................................................................................... 112 Fig. 3.45 2011 comparison of final balanced water surface elevation to the observed elevation at Bankhead Dam.................................................................................... 113 Fig. 3.46 2011 comparison of modeled and observed water surface elevations at Cordova USGS station ......................................................................................................... 114 Fig. 3.47 2011 comparison of final balanced water surface elevation to the observed elevation at Cordova USGS station ........................................................................ 115 xiii List of Abbreviations ANSWERS Areal Nonpoint Source Watershed Environmental Response Simulation APC Alabama Power Company BASINS Better Assessment Science Integrating Point and Non-point Sources BMP Best Management Practice CSTR Continuously Stirred Tank Reactor CWA Clean Water Act CWIS Cooling Water Intake Structure DEM Digital Elevation Model EFDC Environmental Fluid Dynamics Code EPA Environmental Protection Agency EPRI Electric Power Research Institute HSPF Hydrologic Simulation Program - Fortran HUC Hydrologic Unit Code ILWAS Integrated Lake Watershed Acidification Study LAI Leaf Area Index LNSE Logarithmic Nash-Sutcliffe Coefficient NED National Elevation Dataset NLCD National Land Cover Dataset NHD National Hydrography Dataset xiv NRCS National Resources Conservation Service NSE Nash-Sutcliffe?s Efficiency PBIAS Percent Bias RSR Ration of the RMSE to the Standard Deviation SSURGO Soil Survey Geographic Database SWAT Soil and Water Assessment Tool SWMM Storm Water Management Model TMDL Total Maximum Daily Load USGS United States Geological Survey WARMF Watershed Analyses Risk Management Framework WASP Water Quality Analysis Simulation Program 1 Chapter 1 Introduction and Project Purpose 1.1 Background Developing a hydrodynamic model is a time-consuming task that requires extensive field data for set up and calibration. In order to properly represent a system, all of the inputs to that system must be accurately quantified. If poor quality data are input to the model, then a good calibration cannot be obtained or physically inappropriate coefficients must be input to obtain a good calibration. Additionally, each model has built-in assumptions. If the application does not fit the model assumptions, then a physically meaningful representation is not feasible. The model that most closely adheres to the particular project should be used. 1.2 Scope and Objectives A hydrodynamic model was previously developed for the Black Warrior River for the segment between Lewis Smith Dam and Bankhead Dam. This segment of the river is commonly called Bankhead Reservoir. The William H. Gorgas electric generating facility, managed by the Alabama Power Company (APC), is located on this reservoir, approximately halfway between Smith and Bankhead dams. The Gorgas plant is a coal-fired facility equipped with once-through cooling. The hydrodynamic model was developed using EFDC (Environmental Fluid Dynamics Code) (Hamrick 1992). The EFDC model had difficulties predicting overall trends in the reservoir. The inaccuracy of the EFDC model is thought to be a result of a poorly quantified mass balance of the system. There are multiple tributaries which flow into the Black Warrior 2 River over the 75-mile domain which are not accounted for in the model. There is also a known leakage issue at Bankhead Dam. Therefore, the purpose of this thesis is to account for all of the inflows and outflows of the system and evaluate whether adding the resulting mass balance to the EFDC model results in a more accurate model. Certain inflows and outflows from the system are known, such as the inflow from Smith Dam, the discharge (non-leakage) through Bankhead Dam, the withdrawal intake for Plant Miller, and the drinking water intake for the Birmingham Water Works. Also, the four largest tributaries have flow data measured at USGS gage stations, but these stations are located significantly upstream from the reservoir, therefore their data do not represent the true flow entering the domain. Other inflows are not easily quantified and require estimation or modeling. There is an unsteady, unknown leakage rate at Bankhead Dam as well as the water loss due to locking. Additionally, the primary unknown flow contribution is the watershed runoff, which includes the tributary inflows. In order to represent these flows, a watershed model was developed using the watershed model WARMF (Watershed Analysis Risk Management Framework) (Chen et al. 2001). The watershed model is driven by the hydrologic model, which models the movement of water within the watershed as a response to atmospheric input. The hydrologic modeling approach taken with WARMF differed from the traditional hydrologic modeling approach. Typically, hydrologic models are developed to predict peak flows and also to serve as predictive tools. In this application, the time period of interest was over the summer dry season of two years, 2010 and 2011. A traditional split-sample calibration and validation procedure was first followed in an attempt to generate a general 3 hydrologic model (Chapter 2). However, this model will be shown to be inadequate for the desired output of low-flow summer streamflows. Therefore, two separate models presented in Chapter 3 were developed for the summer periods of 2010 and 2011. These highly optimized hindsight (hindcast) models were targeted for the specific purpose of predicting the watershed streamflows for a particular time domain. Therefore, they are not applicable for either another time domain or for predictive modeling. This approach is not traditional, but achieved the goals of modeling streamflows in the watershed for specific time domains. Evaluating the water balance of the system was accomplished using the water balance utility included as part of the 2-D hydrodynamic model CE-QUAL-W2 (Cole 2011). The hydrodynamic model itself was not used. The water balance utility is a stand-alone tool designed to be coupled with the hydrodynamic model and use the same model-control input files for CE-QUAL-W2. It requires a reservoir bathymetry data file developed for CE- QUAL-W2. In this application, the utility was manually coupled to EFDC to determine a mass balance for the reservoir. The utility compares the observed water surface elevation to the original modeled water surface elevation (from either CE-QUAL-W2 or EFDC). The utility then outputs a time-series flow file that when included in the hydrodynamic model will theoretically make that modeled water surface elevations match the observed elevations. If the modeled water surface elevations match the observed elevations, then a mass balance of the reservoir is achieved. The objectives of the study can be outlined as follows: 1) Development of a traditional calibrated and validated hydrologic model for the Bankhead watershed 4 2) Development of targeted hydrologic models to predict streamflows for the specific time periods in 2010 and 2011 3) Development of water balance for the reservoir using the watershed input and the CE- QUAL-W2 water balance utility and application of the balance to the EFDC model to evaluate model improvement 1.3 Study Area The Bankhead watershed is the area that drains into Bankhead Reservoir, between Smith and Bankhead Dams on the Black Warrior River in northwest Alabama. The length of the river is approximately 75 miles from Smith to Bankhead Dam. Figure 1.1 gives a spatial reference where the reservoir is located. The urban center of Birmingham, AL is located in the center of Jefferson County, with the smaller cities of Jasper and Cullman located in central Walker and Cullman counties, respectively. Figure 1.2 shows the location of towns where meteorologic data were available. Figure 1.3 shows the location of the tributaries flowing into the Black Warrior River and Bankhead Reservoir. Figure 1.4 gives the locations of the four USGS gage stations located in the watershed that can be used for model calibration, in addition to the USGS gage on the Black Warrior River (Bankhead Reservoir). 5 Fig. 1.1 Overview of Bankhead Reservoir in relation to neighboring counties Fig. 1.2 Location of available meteorological data in towns near Bankhead Reservoir 6 Fig. 1.3 Tributaries contributing to the Black Warrior River and Bankhead Dam Fig. 1.4 Locations of USGS gage stations in the Bankhead watershed 7 1.4 Thesis Organization Chapter 1 is an introduction to the purpose and goal of the thesis. Chapter 2 gives background information on the WARMF model, details the process of setting up the model, and explores the calibration and validation procedures for the development of the traditional hydrologic model. Chapter 3 outlines the setup of the specialized low-flow WARMF calibrations, the setup of the CE-QUAL water balance utility, and compares the EFDC models with and without the water balance. Chapter 4 is a summary of each step of the process and a recommendation for future study. 8 Chapter 2 Development of WARMF Model 2.1 Introduction WARMF (Chen et al. 2001) is a physically based, versatile, dynamic watershed model that simulates hydrology and water quality. Originally developed for EPRI (Electric Power Research Institute) by Systech Engineering, WARMF has been endorsed by the EPA as a watershed assessment model. WARMF has been compared to other watershed models (Chen and Herr 2002; Chen et al. 2005; Neilson et al. 2003), passed multiple peer reviews (Keller 2000; Keller 2001), used for TMDL development (Chen et al. 2000; Herr et al. 2002; Herr et al. 2003) and nutrient management strategy (NC DENR 2009), and has also been applied for study of future climate scenarios (Rich et al. 2005; Shrestha 2010). WARMF is an integrated watershed model that groups simulation models and input databases into a GIS- based user interface. WARMF was developed from well-established codes, such as ILWAS (Chen et al. 1983), ANSWERS (Beasley et al. 1980), SWMM (Huber et al. 1988), and WASP (Ambrose et al. 1993). WARMF was chosen for this project because of its ability to simulate point and non- point sources, integration of stream and reservoir models into a seamless river basin model, and ability to incorporate spatially varying land uses and meteorological data. WARMF also provides the flexibility to run on an hourly or daily time step, as opposed to other common models such as SWAT (Soil and Water Assessment Tool) (Arnold and Soil 1994) which run on daily time steps. Additionally, WARMF is not constrained by HUC units or an EPA reach 9 file as is HSPF (Chen and Herr 2002), but can delineate the watershed of interest from a DEM. For this study, a WARMF model was set up for the Bankhead Reservoir Watershed to examine the tributary inflow into the Bankhead Reservoir. Because the primary focus of this study is the hydrologic function of WARMF, no water quality study was attempted. However, once the hydrologic model is calibrated, then the WARMF model can be implemented further if water quality studies are desired. In order to set up the model, the watershed was delineated using BASINS (USEPA 2004), input data were compiled from multiple sources, and the hydrologic parameters were calibrated and validated for flow. 2.2 Model Description WARMF is an integrated watershed model and decisions support system that groups models, databases, and decision support tools under one GIS-based framework. It is currently available in the public domain from the U.S. EPA website. WARMF functions through five linked modules: Engineering, Data, Knowledge, Consensus, and TMDL. The modular system is designed as a hub for multiple stakeholders to evaluate different scenarios and reach a decision regarding a Best Management Practice (BMP). For this study, the modules used were the Engineering module, which performs the modeling work, and the Data module, which serves as the database for input files into the model. The Engineering Module of WARMF contains a dynamic watershed simulation model that calculates daily or hourly surface runoff, groundwater flow, non-point source loads, hydrology, and water quality of river segments and stratified reservoirs. The watershed of interest is divided into smaller catchments to better represent the spatial variation within 10 the watershed. Within each catchment, coefficients and physical parameters are constant. Therefore, WARMF lies between being a distributed model and a lumped model. WARMF divides the watershed of interest into a seamless, connected series of catchments, streams, and reservoirs. The publicly available version of WARMF allows for up to 100 catchments. Catchments are further divided into layers, shown in Figure 2.1. Each catchment can have up to five soil layers characterized by their thickness, hydraulic conductivity, saturated moisture content, and field capacity. Land use is aggregated by catchment, which determines the amount of canopy and impervious soils in each catchment. Each layer is modeled as a continuously stirred tank reactor (CSTR) (Chapra 1997). The conceptual hydrologic model within WARMF is shown in Fig. 2.1. The model Fig. 2.1 Hydrologic components within WARMF (Chen et al. 2001) simulates each layer as a hydrologic compartment. The amount of canopy interception is determined by Equation 2.1, which varies as the leaf area index (LAI) changes from month to month. I(mon) is the potential canopy interception for a month, Dmax is the maximum 11 interception possible, L(mon) is the leaf area index for that month, and Lmax is the maximum leaf area index. ( ) ( ) (2.1) After canopy interception is determined to be at the maximum, the excess reaches the snowpack, and eventually reaches the top soil layer. After the available depression storage is filled, the excess water is treated as a sheet flow to the nearest river segment if the soil is frozen or saturated. Sheet flow is predicted by Manning?s equation (Equation 2.2), where Qs is runoff from pervious soils, Z0 is the water depth available for sheet flow, n is the Manning?s roughness coefficient, S is the topographic slope, and where W is the width of the catchment parallel to the receiving stream (Chen et al. 2001). (2.2) If the soil is not saturated or frozen, then the moisture is absorbed into the top layer of soil. From there, the groundwater flows according to Darcy?s Law (Equation 2.3), where Qj is lateral exfiltration from layer j, Khj is the hydraulic conductivity of layer j, and Zj is the thickness of layer j (Chen et al. 2001). (2.3) Each layer from the canopy through the snowpack is treated as a CSTR for mass balance (Equation 2.4) and flow routing (Gherini et al. 1985). For each time step, the model performs a water balance to update soil moisture contents after evaporation, exfiltration, and percolation between soil layers. For each soil layer, the water balance is given in Equation 2.3, where Vj is the volume of water in the soil layer j, Vj0 is the initial volume of water in the soil layer j, Ij is the infiltration to layer j. Ij+1 is the percolation from layer j to layer j+1, Lj is 12 the lateral inflow from an upstream segment, Ej is the evapotranspiration from layer j, and Qj is the lateral exfiltration from layer j (Chen et al. 2001). (2.4) If Vj from Equation 2.4 is greater than the volume of the saturated moisture content, then the infiltration is impeded, and the excess is returned to the layer above. Thus the model computes mass balance from the bottom soil layer and moves upward. Stream segments are routed by the kinematic wave method from the uppermost tributaries to the terminus of the watershed. For each segment, the change in storage is the inflow minus the outflow (Equation 2.5). The outflow from the streams is determined by Manning?s Equation, shown in Equation 2.6 (Chen et al. 2001), where O is the outflow, D is the hydraulic radius of the stream, S is the bed slope, A is the cross-sectional area of the stream, and n is Manning?s roughness coefficient. (2.5) (2.6) The potential evaporation is the maximum evaporation from the free water surface and the transpiration from the soil that can occur. If insufficient water is available, then only the available water will be transpired. The potential evaporation is calculated for each month according to latitude by the Hargreaves formula(Gao et al. 2008). The monthly potential is converted to a daily potential by Equation 2.7. ( ) ( ) (2.7) For water quality simulations, the model simulates wet deposition, dry deposition, fertilizer, pesticide, animal droppings, dissolution, advection, decay, soil erosion and wash 13 off, mineral weathering, acid mine drainage, septic systems, organic matter decay, nitrification, cation exchange, and plant uptake. 2.3 Literature Review The hydrologic cycle is comprised of six major components: precipitation, evaporation, transpiration, infiltration, surface runoff, and groundwater flow ( Fig. 2.2). The governing physical principles behind the cycle are constant, but different modeling methods approach each element differently. Pollutants (such as nutrients) can also be modeled as part of the cycle. The movement of water from component to component affects the pollutant behavior and how it is transported. Surface runoff will carry pollutants, whereas evaporation will leave pollutants behind. TMDL management arose out of a need to establish Best Management Practices (BMPs) with regard to nutrients and other pollutants. All sources of water and pollution can be placed into two categories: point-source discharges, and non-point-source discharges. Point discharges are defined by the discharge of water or a pollutant at a single, finite location, such as a pipe. Non-point-source discharges are spread out over an area typically as an atmospheric deposition, of which the primary case is rainfall. The surface runoff from rainfall picks up pollutants, which are deposited into surface waters. When choosing a hydrologic model, it is important to choose a model that meets the particular needs of a project (Parsons et al. 2004). There are dozens of models available, but each model tends to focus on a particular intended use. A few of the more commonly used watershed models are described below. 14 Fig. 2.2 Representation of hydrologic cycle (Dardashti 2010) 2.3.1 BASINS/HSPF (Hydrologic Simulation Program ? Fortran) BASINS (Better Assessment Science Integrating Point and Non-point Sources) is a model developed by the EPA for determining TMDLs. The hydrologic component within BASINS is HSPF, Hydrologic Simulation Program ? Fortran (Donigian et al. 1984). Like WARMF, BASINS/HSPF is applicable to large watersheds, however, BASINS is inextricably tied to the resolution of HUC (Hydrologic Unit Code) and the EPA Reach File (Chen and Herr 2002). WARMF is capable of being used on a watershed of any size, because the only required spatial input is the Digital Elevation Model (DEM). HSPF, similar to WARMF, is able to run on an hourly time step (USEPA 2004). WARMF is a mechanistic model, while HSPF is an empirical water budget model (Chen et al. 2005). HSPF has a temperature adjustment factor that translates the temperature at the monitoring station to the 15 subwatershed. However, it does not have an adjustment factor for precipitation (Chen and Herr 2002) as WARMF does. In a direct comparison between BASINS and WARMF on the Mica Creek Watershed, HSPF was forced to get rid of excess water by sending 30-40% to an inactive groundwater sink. WARMF did not have a similar issue, due to its ability to weight precipitation factors (Chen et al. 2005). While BASINS is designed to be linked to HSPF, it is a tool useful for other applications. Built on the open-source MapWindow GIS platform, it contains tools for various spatial analyses. It is capable of manual and automatic watershed delineations, which can then be exported into other programs. 2.3.2 SWAT (Soil and Water Assessment Tool) The Soil and Water Assessment Tool (SWAT), was developed for the USDA (United States Department of Agriculture) with the goal of evaluating the impact of various management practices on large watersheds (Arnold and Soil 1994). Unlike the mechanistic hydrologic model in WARMF, SWAT employs the NRCS curve number method to determine peak runoff rates (Parsons et al. 2004). The method for computing subsurface flow is a kinematic storage routine (Neitsch et al. 2005) that does not consider detailed information of a tile-drain system. SWAT can be run on a monthly or daily time step. SWAT has been implemented successfully (Bingner 1996), but it has also shown a tendency to underestimate flows in the spring and overestimate flows in the fall (Gollamudi 2006). Additionally, the model has struggled to predict exceptionally wet years (Chu and Shimohammadi 2004). SWAT has also been shown to give less accurate results for shorter time intervals (daily) rather than longer (monthly) time intervals (Dardashti 2010). 16 Both SWAT and WARMF are physically-based, distributed-parameter and continuous watershed modeling tools. An application of both WARMF and SWAT for TMDL management (Wang et al. 2004) revealed the following limitations of each model in regards to each other. WARMF requires more computational effort than SWAT, which is a result of the complex governing equations used in WARMF along with the higher spatial resolution. WARMF tracks minerals which affect the fate and transport of nutrients that SWAT does not model. However, WARMF does not clearly separate lateral flow and ground water return flow, and does not make simulated values for a number of parameters available to users. Unlike WARMF, SWAT does not include an input parameter for septic systems as a nutrient source nor does it allow for dynamic air quality input data as a source of pollutants. Also, SWAT was developed for agricultural purposes and is not suitable for modeling urban land uses. An advantage of SWAT is its ability to accept GIS-based soil data, rather than requiring soil data to be input manually (Geza and McCray 2008). The accuracy of the hydrological model within SWAT has been sufficient for its intended purpose of developing long-term land-use management strategies, but when directly compared to WARMF, SWAT predicted the observed values less accurately. For a study conducted on the Napa River, the Nash-Sutcliffe values were 0.59 and 0.75 for SWAT and WARMF, respectively (Wang et al. 2004). 2.3.3 WARMF WARMF was originally developed under sponsorship from EPRI, by Systech Engineering as a decision support system for TMDL management of watersheds. The model has been successfully applied to numerous watersheds throughout North America and the 17 world (Luo et al. 2011). WARMF has undergone extensive testing, and been compared with similar watershed models in multiple studies. WARMF was evaluated by a panel of independent experts according to the ?Guidelines for Conducting External Peer Review of Environmental Regulatory Models? (EPA 100-B-94-001). At the time of its review, WARMF was the only watershed-scale model that had been so extensively reviewed (Keller 2000). Algorithms for different components within WARMF are taken from established codes. The primary computational base is from the Integrated Lake-Watershed Acidification Study (ILWAS) model (Chen et al. 1983), which divides a watershed into catchments, stream segments, and reservoir layers (Chen et al. 2001). Flow is routed through the canopy, snow pack, into the soil layers, exfiltration from the soil, evaporation from the soil, through stream segments by the kinematic wave method, through soil layers by Darcy?s Law, and through stratified reservoirs to the terminal point of the watershed. The original ILWAS model was only capable of modeling a single terminal reservoir, but WARMF was modified to model multiple reservoirs. The algorithms for modeling pollutant transfer and soil erosion from land surfaces were taken from ANSWERS (Beasley et al. 1980). The algorithm for modeling pollutant wash-off from urban land uses was adapted from SWMM (Huber et al. 1988). The kinetics of nutrients, algal dynamics, and sediment sorption-desorption of pesticides and phosphorus were adapted from WASP (Ambrose et al. 1993). WARMF has been used to model various scenarios with success. Various studies include modeling TMDLs in the Chartiers Creek watershed (Chen et al. 2000) and Hangman Creek watershed (USEPA 2007), modeling onsite wastewater systems in Turkey Creek 18 Watershed (Geza 2009), modeling DO in the San Joaquin River (Stringfellow et al. 2009), Falls Lake watershed (NC DENR 2009), and modeling acid mine drainage in the Cheat River Basin (Chen et al. 2004). Initial applications of WARMF were in cooperation with EPRI members and included the Catawba River Basin, Holston River Basin, Hockanum River Basin, Cheat River, Oostanaula Creek, Dillon Lake Watershed, and Mica Creek (Chen et al. 2001). A particular strength of WARMF is its flexibility to be applied to different scenarios and even modified to extend its functionality beyond what was originally intended. WARMF has been used to model the effect of septic outflow (Geza and McCray 2007; Weintraub et al. 2004), applied to future climate change scenarios (Dayyani et al. 2012; Shrestha 2010), used to analyze the difference between observed and modeled atmospheric data (Herr et al. 2010), modified to support the ZeroNet water-energy initiative (Riggs et al. 2005), modified to work alongside DRAINMOD (Dardashti 2010), and modified to model mercury fate and transport (Chen 2006). WARMF was chosen for this project because it fit the needs of the project more closely than other models. The primary goal of this study was to achieve a good prediction for flow, and WARMF has produced better calibrations for flow rates than some other models (Wang et al. 2004). Also, the Bankhead watershed includes urban areas, which SWAT is not appropriate to model. The watershed includes areas from several HUC areas, which is difficult for some models to incorporate. Finally, extensive data were available for input, which made an hourly time step in WARMF feasible. The most significant limitation of WARMF is that the publicly available model is limited to 100 catchments. For hybrid lumped and distributed models, this limitation can be 19 significant for large watersheds, because larger areas must be lumped together and the physical parameters across that larger land area are averaged out. For example, a mountainous area and a more plateau-like area could be averaged together, rather than treated as separate catchments. Smaller catchments result in more physically meaningful discretization of the watershed 2.4 Methods 2.4.1 Watershed Delineation Using Basins The principle required input for WARMF is a delineated watershed that includes land catchments, stream segments, and reservoirs (if present). WARMF has the ability to delineate a watershed from an input DEM (Digital Elevation Model), but this feature is disabled in the publicly-available version of WARMF. It is recommended to use EPA?s BASINS software to create the watershed, which can be then imported into WARMF (Systech Engineering Inc 2005). BASINS is an EPA-distributed software based on a GIS platform that has built-in connections to many national datasets. It incorporates an open-source GIS system, MapWindow, as the GIS engine underneath the data. The download tool within BASINS downloads data such as 8-digit HUC units, NHD (National Hydrography Dataset), DEM, NED (National Elevation Dataset), NLCD (National Land Cover Dataset), and various political boundaries (cities, states) and census data. The Bankhead Watershed was determined to be composed of two complete HUC units (Mulberry Fork and Locust Fork) and two partial HUC units (Sipsey Fork and Upper Black Warrior), shown in Fig. 2.3. 20 In order to create the delineation shown in Fig. 2.3, a DEM covering the area with 30- meter resolution was downloaded via the download tool in BASINS. A river network can be downloaded to serve as a guide for the delineation process, but was found to be unnecessary, as the delineations with and without the river network were identical. Data from the NED (National Elevation Dataset), with a resolution of 10 m, are also available, but are not recommended for use with watersheds of this size (Systech Engineering Inc 2005). The finer resolution can potentially have many small ?holes? in the terrain which do not drain to the rest of the watershed. These holes are artificially filled during the delineation, but can still cause issues. Fig. 2.3 HUC units contributing to Bankhead Watershed 21 The BASINS delineation tool allows for certain points to be chosen that will become the downstream point of a catchment, which is useful if USGS stream gage data are available. For the delineation, the four available USGS stations were set to be outflow points for catchments in order to compare modeled and observed data from the same location. The completed delineation resulted in 83 catchments (the publicly available version of WARMF allows a maximum of 100 catchments), with a total area of 3,064 mi2, or 7,849 km2, with an average area of 36.9 mi2, or 94.6 km2. The slopes of the catchments ranged from 1.1% to 10.1%, with a mean of 4.7% (Fig. 2.4). Fig. 2.4 Average slopes of delineated catchments In order to incorporate the reservoir into the catchment system, the river segments associated with the catchments on the reservoir were deleted and replaced with reservoir 22 segments. Together, the land catchments, river segments, and reservoir are the layers input into WARMF. WARMF accepts land cover data as a shapefile input, so the NLCD 2001 was downloaded from BASINS, converted to a compatible shapefile format, and input into WARMF, which automatically calculates the percentages of land cover types represented in the watershed. A summary of the land cover found in four sample segments is shown in Table 2.1. The locations of these sample catchments are shown in Fig. 2.6.. Fig. 2.5 Layout of delineated catchments with river segments and reservoir in WARMF 23 Table 2.1 Land Cover data from sample catchments Land Cover Subcatchment 3 Subcatchment 76 Subcatchment 57 Subcatchment 11 (Agricultural region) (Birmingham metro area) (Gorgas plant area) (Forested region) Percentage Percentage Percentage Percentage Deciduous 18.95 14.53 33.8 22.44 Coniferous 3.62 8.54 29.86 25.88 Mixed Forest 4.98 2.88 9.69 10.18 Orchard 0 0 0 0 Cropland / Pasture 63.93 4.74 12.24 27.11 Confined Feeding 0 0 0 0 Rangeland 0.95 0.94 4.77 6.27 Forested Wetland 0.39 0.85 1.03 2.13 Non-forested Wetland 0.01 0 0.12 0.02 Tundra 0 0 0 0 Barren 0.05 0.28 2.1 0.35 Residential 6.33 49.66 2.47 5.09 Comm./Industrial 0.57 17.05 0.05 0.09 Water 0.22 0.53 3.87 0.42 TOTAL 100 100 100 100 24 Fig. 2.6 Location of sample catchments for land cover data 2.4.2 WARMF Model Input WARMF requires input parameters of time-series data to simulate watershed hydrology and water quality. The required inputs are either hourly or daily time series of atmospheric data and air quality data. The air quality data are not required to be on a continuous, regular time step and is only important if water quality is to be simulated. For this project, since water quality simulation was unimportant, a ?blank? air quality file was used. 25 2.4.2.1 Observed Streamflow Data Observed streamflow time series are necessary for model calibration. In the watershed, there are four USGS gage stations (Table 2.2) with data available for download. There is an additional gage station located on the Black Warrior River near Cordova, but this section of river is set up as a reservoir, so it was not used in the calibration of the WARMF model, but was used later to calibrate the hydrodynamic model. Table 2.2 USGS gage stations located in Bankhead watershed Gaging Station No. Stream Data Record Period USGS 02453000 Blackwater Creek near Manchester, AL (drainage area 487 mile2) Oct 1, 1938 to present USGS 02450180 Mulberry Fork Near Arkadelphia, AL (drainage area 181 mile2) Oct 1, 1976 to present USGS 02456500 Locust Fork at Sayre, AL (drainage area 258 mile2) Oct 1, 1928 to present USGS 02454055 Lost Creek above Parrish, AL (drainage area 143 mile2) Oct 1, 1992 to present 2.4.2.2 Atmospheric Data The required meteorological input data for WARMF are precipitation (cm), minimum temperature (?C), maximum temperature (?C), cloud cover, dew point temperature (?C), air pressure, and wind speed. Hourly data were available from the Alabama Mesonet Data program, maintained by AWIS. The closest available station was located in Cullman, AL, in the northeast corner of the watershed. Some data required processing. Cloud cover is not provided by AWIS, so it was calculated by the method in Equation 2.7 (Systech Engineering Inc 2005), where Tmin is the minimum temperature, Tmax is the maximum temperature, Tdew is the dew point temperature, Precip is the precipitation, and Cloud is the resulting cloud cover. 26 The algorithm estimates the cloud cover based upon the relationship between the temperature and the dewpoint temperature at each timestep if there is no precipitation. Smaller differences between air temperature and dewpoint temperature result in an estimation of higher cloud cover. If precipitation occurs at a timestep, then the cloud cover ranges from 0.8 to 1 (full cloud cover) based upon the intensity of the precipitation. (2.7) Dew point temperature was also not directly measured, but calculated from the measured temperature and relative humidity values. Dew point temperature was calculated according to the following method in Equation 2.8 (Buck 1981), where RH is the measured relative humidity, Tdp is the dew point temperature, and b and c are constants equal to 17.368 and 238.88, respectively. ( ) ( ) (2.8) 27 APC provided rainfall data from twelve separate sites scattered throughout the basin area, shown in Fig. 2.7. These sites only provided rainfall on a daily time step, when hourly data were needed as input. In order to dis-aggregate the daily data to an hourly time step, the hourly precipitation data from Cullman Mesonet were converted to a percentage of daily rainfall that occurred in each hour. These daily distributions were then applied to the daily rainfall totals measured at the other locations to get hourly rainfall data. Therefore, the total daily precipitation for each day was left unchanged, but spread out over the day according to a distribution based on the pattern observed at the Cullman Mesonet site for that day. Fig. 2.7 Location of meteorological stations in watershed region 28 A meteorological time series must be defined for each catchment. Catchments were assigned to the closest available meteorological station. 2.4.2.3 Soil Data Soil parameters are a key input for each catchment. Soil data were downloaded from the SSURGO database maintained by the NRCS. Soil data are provided by county. Select soils from each county are shown in 2.4.3 WARMF Simulation The simulation for WARMF was set to run on an hourly time step from January 2, 2010, to September 1, 2010. This time period was used as the calibration period to tweak the input parameters. After achieving the best calibration, the time period from September 2, 2010 to June 30, 2011 was used to validate the model calibration. This time frame allowed the model to be evaluated under a variety of conditions, including periods of heavy precipitation as well as the seasonal dry period in late summer. Table 2.3. The soils shown are the most prevalent soil types in the county. Each county has numerous soil types present in small percentages. The SSURGO database gives values for the soil layer thicknesses as well as the maximum hydraulic conductivity. These values serve as a starting point for the WARMF calibration, but will be tweaked as part of the calibration process. Table 2.3 shows how soil parameters can vary widely from one soil type to another. The process of lumping catchments together groups these disparate soils into one homogenously parameterized soil. As catchment sizes increase, more soils are lumped into this parameter, and more physical reality is lost. 29 2.4.3 WARMF Simulation The simulation for WARMF was set to run on an hourly time step from January 2, 2010, to September 1, 2010. This time period was used as the calibration period to tweak the input parameters. After achieving the best calibration, the time period from September 2, 2010 to June 30, 2011 was used to validate the model calibration. This time frame allowed the model to be evaluated under a variety of conditions, including periods of heavy precipitation as well as the seasonal dry period in late summer. Table 2.3 Soil data downloaded from SSURGO for the three major counties 30 2.4.4 Calibration Criteria Evaluating the accuracy of a hydrologic model includes two key approaches: graphical inspection and statistical evaluation. Graphical evaluation compares the general shape of the plots of observed values and modeled values to evaluate the ?fit? of the model. This technique is useful for determining whether the peak flows are modeled accurately and 31 whether the modeled values follow the trends of the observed values. However, graphical techniques are unable to give a quantitative measure of the fit of a model and determine which parameters give a more statistically accurate prediction. 2.4.4.1 Literature Review of Hydrologic Model Calibration The most common traditional measures of the accuracy of a model are the coefficient of determination (R2) and the RMSE (Root Mean Squared Error) value. Shortcomings of these measures led to the development of new measures, such as the Nash-Sutcliffe Efficiency (NSE), ratio of the RMSE to the standard deviation (RSR), and the Percent Bias Ratio (PBIAS). Built-in mathematical biases in the equations for these parameters mean that there can be no perfect statistic that evaluates the model without bias and gives the best possible ?fit? under every situation. However, using multiple parameters to evaluate the model rather than maximizing one particular parameter should best reflect the dynamics of the entire system (Krause et al. 2005). Equations showing the method for calculating the previously mentioned parameters are shown below in Equations 2.9-2.12, where O are observed values, P are predicted values, and i represents the i-th time step. 2.4.4.2 Coefficient of Determination (R2) The coefficient of determination (R2) is defined as the squared value of the coefficient of correlation, or as the squared ration between the covariance and the multiplied standard deviations of the observed and predicted values (Krause et al. 2005). R2 is calculated according to Equation 2.9 and ranges from 0.0 to 1.0, with higher numbers indicating better fit. R2 should not be used as the sole measure of a model?s accuracy because it is insensitive 32 to the proportional differences between the simulated and observed values and can still give high values even if the modeled and simulated values are significantly different in magnitude (Legates and McCabe 1999). Additionally, R2 is highly sensitive to outlying values that greatly differ from the mean, thus a model which accurately predicts peak events but poorly predicts values nearer to the mean will have an inflated R2. ? ( ?)( ?) ?? ( ?) ?? ( ?) (2.9) 2.4.4.3 Nash-Sutcliffe Efficiency (NSE) The Nash-Sutcliffe coefficient (E) is defined as one minus the sum of the absolute squared differences between the predicted and observed values normalized by the variance between the observed value and the mean observed value (Krause et al. 2005). It is computed using Equation 2.10 (Nash and Sutcliffe 1970). ? ( ) ? ( ?) (2.10) The value for E ranges from -? to 1, with higher values indicating better agreement (Legates and McCabe 1999). Values below 0 indicate that the average observed value is a better predictor than the model. The Nash-Sutcliffe coefficient is an improvement over R2 because it is more sensitive to the difference between the modeled and observed values. However, it is still highly sensitive to outliers due to the squared terms. Despite its similarity to R2, E has been shown to be un-correlated to R2 (Krause et al. 2005). 33 Because E has become widely used for evaluating hydrologic models, attempts have been made to improve upon its formulation to reduce its sensitivity to extreme values and increase the relative influence of low flows (Krause et al. 2005). The most straightforward modification is taking the natural logarithm of the observed and predicted values before calculating E. This transformation has the effect of reducing peak flows, but not heavily modifying low flows, therefore increasing the influence of low flows at the expense of the influence of extreme points. A second modification involves using the absolute difference of the observed and predicated values rather than squaring the difference. This transformation has the same effect as the logarithmic approach, because the huge numbers resulting from the squaring of large numbers are not involved. The resulting modified Eabs is always less than the unmodified E, which allows for a larger calibration range, but can give the impression that the model is less accurate than it actually is (Krause et al. 2005). 2.4.4.4 RMSE and Standard Deviation Ratio (RSR) RMSE (Root Mean Squared Error) is an established statistic to measure model fit, where the lower value indicates better model performance (Chu and Shimohammadi 2004). The error index provided by the RMSE can be normalized by dividing it by a normalization term, which is the standard deviation of the observed values. Thus the resulting statistic, RSR, can be applied to various constituents (Moriasi et al. 2007). RSR ranges from 0, perfect fit, to a large value. RSR is calculated according to Equation 2.11. ?? ( ) ?? ( ?) (2.11) 34 2.4.4.5 Percent Bias (PBIAS) The Percent Bias estimator differs from the other parameters in that it provides insight as to whether the model has a bias toward overpredicting or underpredicting the flows. PBIAS is calculated according to Equation 2.12, and the ideal value is 0.0. Positive values indicate a bias toward underestimation and negative values indicate a bias toward overestimation (Gupta et al. 1999). An inherent weakness in the PBIAS statistic is that high error in one portion of the model reflects very poorly in the final statistic. For example, an otherwise perfect model with one section of significant overprediction due to inaccurate data input would result in a large, negative PBIAS value. This high PBIAS leads to the conclusion that the model is very poor, while in reality, the model is only inaccurate in that one time domain. ? ( ) ? *100 (2.12) 2.4.4.6 Evaluation Statistics Used For this project, the suggested three parameters of Moriasi et al. (2007) of Nash- Sutcliffe Efficiency (E), RSR, and Percent Bias were used to evaluate the hydrologic model developed in WARMF. Nash-Sutcliffe Efficiency (unmodified) was selected for its consistent use in the literature, which provides an established benchmark of model performance. The RSR value was selected because it includes an error component (RMSE) as well as the normalization benefit of the standard deviation (Moriasi et al. 2007). PBIAS was also used to expose tendencies in the model to overpredict or underpredict streamflows. 35 Graphical and statistical analyses must be used simultaneously to evaluate a model. A time-series plot of discharge vs. time is the most commonly used graph in evaluating a hydrologic model. This approach is primarily useful for visualizing how accurately the model predicts peak flows. Also, particular sections of the graph can be analyzed to determine the predictive accuracy of baseflow, as well as the rising and fall limbs of the hydrographs for each storm event within the time period. However, this approach alone is not sufficient. A model may appear to be a very good fit from this perspective, but does not easily show consistent over or under estimation of flows. A graph which plots the cumulative flow over the simulated time period will reveal consistent bias within the model. Also, this type of graph can be used to analyze how well the model predicts base flow by comparing the slope of the cumulative lines, as well as peak events by comparing the relative jumps in the cumulative graphs. The graphical techniques discussed are necessary to evaluate a model?s performance, along with the statistical parameters already discussed. Blindly optimizing a single statistical measure does not necessarily optimize the model. For example, a PBIAS of 0 could be obtained by using the mean flow for the time period. This approach is clearly not a valid model, but that statistic used blindly would suggest that it is a perfect model. Also, graphical analyses typically reveal the weaknesses in the model and guide changes to improve each area of model?s accuracy. An additional reason that graphical and statistical analyses must be used simultaneously is that graphical comparisons can show that inaccuracies may not be the fault of the model. For example, the observed data may show that a storm event occurred, but the 36 model shows no evidence of rainfall. This particular case might not be a fault of the model, but could be a shortcoming of the input data. Statistical benchmarks for hydrologic modeling (Moriasi et al. 2007) for streamflow are given in Table 2.4. However, the guidelines listed in Table 2.4 are for models that use daily time steps. In model evaluation, shorter time steps tend to lead to poorer statistical results (Moriasi et al. 2007). Therefore, the numerical benchmarks employed to evaluate this hourly model could be considered stringent rather than lenient. An additional consideration in model evaluation is the accuracy of the observed streamflow data. USGS gage streamflow has a variable, unknown accuracy for each station, so even a ?perfect? calibration may not reflect the real-world situation perfectly. Table 2.4 Statistical benchmarks for evaluating hydrologic models (Moriasi et al. 2007) Performance Rating NSE RSR PBIAS Very Good 0.75 < NSE < 1.00 0 < 0.50 PBIAS < ?10 Good 0.65 < NSE < 0.75 0.50 < RSR < 0.60 ?10 < PBIAS 0.70 PBIAS >?25 2.5 Results and Discussion The WARMF model was run for two separate time periods on an hourly time step. The calibration period was January 20, 2010 through September 1, 2010. The model was given a spin-up time from January 1 to January 20 before the evaluation period began. After achieving a good calibration, the model was verified using data from September 2, 2010 to June 30, 2011. 37 Each of the four gaged subwatersheds was calibrated separately to produce the best results. The coefficients which yielded the most statistically and graphically accurate results were then analyzed and tweaked so that the coefficients across all the catchments were similar. Some variation among catchments was allowed, based upon the soil type and distribution within the catchments. 2.5.1 Statistical Summary The calibration and validation results were generally ?good? to ?very good,? with the single exception of the PBIAS for the validation of the Blackwater Creek subwatershed, which only graded as ?satisfactory.? The statistical summary of the four watersheds is given in Table 2.5. Below each statistic is a numerical value that represents whether the measure is very good (1), good (2), satisfactory (3), or unsatisfactory (4). Table 2.5 Statistical summary of calibration and validation for streamflow CALIBRATION VALIDATION NSE RSR PBIAS(%) NSE RSR PBIAS(%) Blackwater Creek 0.78 0.47 10.18 0.71 0.54 24.73 1 1 2 2 2 3 Locust Fork 0.88 0.34 -0.01 0.72 0.53 -15.82 1 1 1 2 2 3 Lost Creek 0.79 0.46 -3.00 0.73 0.52 1.91 1 1 1 2 2 1 Mulberry Fork 0.86 0.37 0.16 0.77 0.47 -13.10 1 1 1 1 1 2 The results for each subwatershed are further analyzed and discussed in the following sections. Graphical analysis is incorporated to identify particular strengths and weaknesses of the model. The graphical approaches used are time-series of flow, to examine the model?s 38 predictions of peak flows, as well as cumulative volume, to examine the model?s tendency to have a bias. 2.5.2 Blackwater Creek Blackwater Creek is one of the two smaller watersheds in the model, along with the Lost Creek Watershed. The Blackwater Creek Watershed is composed of the three most northwest catchments, comprising an area of 155 square miles within Walker County, AL. The subwatershed is shown in Fig. 2.8. Fig. 2.8 Blackwater Creek subwatershed 39 The graphical analysis comparing the observed to modeled streamflow (in m3/s or cms) is shown in Fig. 2.9 (calibration) and Fig. 2.10 (validation). The calibration period shows a tendency to underpredict the smaller events, but matches the peak events fairly well. The recession limb of the hydrographs appears to be underpredicted, with the exception of the major events around Julian Day 120 (April 30, 2010). Inspection of the cumulative volume graph confirms that the model has a regular tendency to underpredict the observed values. The trend is followed fairly accurately, and the baseflow during the low-flow period appears to be fairly well-represented. The final volume difference (PBIAS) is 10.2%. The NSE, RSR, and PBIAS are all on the dividing line of being labeled ?good? or ?very good.? As expected, the validation period (Fig. 2.10) does not perform as well as the calibration period. The model continued its trend of underpredicting the observed values. The cumulative graph shows that the large discrepancy in the cumulative volume (PBIAS of 24.7%) is primarily due to significant underprediction near Julian Days 430 and 470 (March 7, 2011-April 16, 2011). The underprediction observed during both the calibration and validation periods could be corrected through adjusting several different model coefficients but primarily by increasing the precipitation weighting factor. However, adjusting the model to optimize the PBIAS results in a lower NSE value. Therefore, the current model coefficients were chosen to give the most balanced model that adequately describes each component of the hydrograph rather than optimizes a single evaluation parameter. The final calibration coefficients are shown in Table 2.6. 40 Table 2.6 Calibration coefficients for the Blackwater Creek subwatershed SOIL LAYERS Layer 1 Layer 2 Layer 3 Prec Weight 0.85 Temp Lapse 3.5 Surface Roughness 0.2 Stream Roughness 0.1 Soil Thickness (cm) 22 25 28 Initial Moisture 0.4 0.4 0.35 Field Capacity 0.35 0.35 0.3 Saturation 0.45 0.4 0.35 Horizontal Conductivity (cm) 190000 150000 100000 Vertical Conductivity (cm) 95000 75000 50000 Root Distribution 0.75 0.1 0 Soil Density 1.3 1.3 1.3 Soil Tortuosity 10 10 10 41 Fig. 2.9 Flow calibration of Blackwater Creek subwatershed 0 20 40 60 80 100 120 140 20 70 120 170 220 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 10000 20000 30000 40000 50000 60000 70000 20 70 120 170 220 Vol um e ( cu bic me ter s) Julian Day Cumulative Observed Cumulative Modeled 42 Fig. 2.10 Flow validation of Blackwater Creek 0 50 100 150 200 250 244 294 344 394 444 494 544 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 10000 20000 30000 40000 50000 60000 70000 80000 244 294 344 394 444 494 544 Vol um e ( cu bic me ter s) Julian Day Observed Flows Modeled Flows 43 2.5.3 Locust Fork Locust Fork is the largest subwatershed in the model. It contains the majority of Blount County and a portion of Jefferson and Cullman counties. The total drainage area is 858 square miles and is the largest tributary flowing into the Black Warrior River. The plots of the calibration and validation periods for Locust Fork are shown in Fig. 2.12 and Fig. 2.13 respectively. Statistically, the model performs extremely well, with an NSE of 0.88 and a PBIAS of -0.01%. Inspection of the calibration plots shows that the model has a tendency to underpredict peak events, but matches the falling limb of the hydrograph well. During the validation period, the model predicted false events for Julian Days 300- 350 (October 28, 2010-December 17, 2010). These predictions can be attributed to precipitation in the input data that the watershed did not actually receive. This influx of flow causes the cumulative curve to become significantly high in this region. From that point, the cumulative curve follows the trends observed values well. The high PBIAS (-15.82%) is primarily due to the error in that timeframe. Overall, the graphical inspections agree with the excellent statistical numbers that the model is very good for this subwatershed. The final calibration coefficients for the Locust Fork subwatershed are given in Table 2.7. 44 Fig. 2.11 Locust Fork subwatershed 45 Fig. 2.12 Flow calibration of Locust Fork subwatershed 0 100 200 300 400 500 600 700 20 70 120 170 220 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 50000 100000 150000 200000 250000 300000 350000 20 70 120 170 220 Vol um e ( cu bic me ter s) Julian Day Observed Modeled 46 Fig. 2.13 Flow validation of Locust Fork subwatershed 0 100 200 300 400 500 600 700 800 244 294 344 394 444 494 544 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 50000 100000 150000 200000 250000 300000 244 294 344 394 444 494 544 Vol um e ( cu bic me ter s) Julian Day Cumulative Observed Cumulative Modeled 47 Table 2.7 Calibration coefficients for the Locust Fork subwatershed SOIL LAYERS Layer 1 Layer 2 Layer 3 Prec Weight 0.9 Temp Lapse 0 Surface Roughness 0.2 Stream Roughness 0.04 Soil Thickness (cm) 15 25 30 Initial Moisture 0.4 0.4 0.35 Field Capacity 0.35 0.3 0.3 Saturation 0.45 0.4 0.35 Horizontal Conductivity (cm) 170000 120000 50000 Vertical Conductivity (cm) 85000 60000 25000 Root Distribution 0.75 0.1 0 Soil Density 1.3 1.3 1.3 Soil Tortuosity 10 10 10 2.5.4 Lost Creek The Lost Creek subwatershed is located directly south of the Blackwater Creek subwatershed, and is also composed of only three catchments. The subwatershed is shown in Fig. 2.14. The total drainage area is 176 square miles. The plots for the calibration and validation periods are shown in Fig. 2.15 and Fig. 2.16, respectively. The calibration period suggests an excellent match of peak flows, with the notable exception that the model completely missed several events. Because the model was an excellent predictor for the other events, these omissions are likely due to input data that failed to represent the actual precipitation. The cumulative curve shows that the model does follow the observed volume well. Interestingly, the missed peak events do not have a huge effect on the cumulative volume, which suggests that cumulative volume is most affected by baseflow conditions. 48 The validation period continues the same trends observed during the calibration period, with a few missed peaks. There is more deviation within the cumulative volume curve, but at the end of the simulation the cumulative volumes become very similar with a resulting PBIAS of 1.91%. The final calibration coefficients are given in Table 2.8. Fig. 2.14 Lost Creek subwatershed 49 Fig. 2.15 Calibration of Lost Creek subwatershed 0 20 40 60 80 100 120 140 160 180 20 70 120 170 220 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 10000 20000 30000 40000 50000 60000 20 70 120 170 220 Vol um e ( cu bic me ter ) Julian Day Cumulative Observed Cumulative Modeled 50 Fig. 2.16 Validation of Lost Creek subwatershed 0 50 100 150 200 250 300 244 294 344 394 444 494 544 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 10000 20000 30000 40000 50000 60000 244 294 344 394 444 494 544 Vol um e ( cu bic me ter ) Julian Day Observed Modeled 51 Table 2.8 Calibration Coefficients for the Lost Creek subwatershed SOIL LAYERS Layer 1 Layer 2 Layer 3 Prec Weight 0.8 Temp Lapse 5 Surface Roughness 0.2 Stream Roughness 0.1 Soil Thickness (cm) 22 25 20 Initial Moisture 0.5 0.5 0.5 Field Capacity 0.35 0.35 0.3 Saturation 0.45 0.4 0.35 Horizontal Conductivity (cm) 45000 45000 45000 Vertical Conductivity (cm) 22500 22500 22500 Root Distribution 0.75 0.1 0 Soil Density 1.3 1.3 1.3 Soil Tortuosity 10 10 10 2.5.5 Mulberry Fork The Mulberry Fork subwatershed is located immediately north of the Locust Fork subwatershed, and is the second-largest tributary of the Black Warrior River. The subwatershed is shown in Fig. 2.17. The drainage area of the Mulberry Fork subwatershed is 491 square miles. The statistics of an NSE of 0.86 and 0.77 suggest that the Mulberry Fork validation and calibration was the most accurate of the four subwatersheds. The calibration period of the Locust Fork subwatershed was marginally higher statistically, but Mulberry Fork was the only subwatershed to have its calibration and validation NSE values qualify as ?very good.? The validation PBIAS was slightly high at -13.1%, which prevented the calibration from achieving a ?very good? rating on every statistic. 52 Graphical analysis in Fig. 2.18 and Fig. 2.19, shows that the model does predict the observed values exceptionally well. The model missed two minor events during the calibration period and underpredicted three minor events during the validation period. However, the cumulative volume curve follows the observed values very well during the calibration period. The cumulative curve for the validation period shows where the model erred and caused the relatively high PBIAS value. The model predicts a small event around Julian Day 340 (December 7, 2010), which is not reflected in the observed values. Thus, this event and the subsequently higher baseflow drive the cumulative curve higher than the observed values. The final calibration coefficients are given in Table 2.9. Fig. 2.17 Mulberry Fork subwatershed 53 Fig. 2.18 Calibration of Mulberry Fork subwatershed 0 100 200 300 400 500 600 700 800 20 70 120 170 220 Flo w (cm s) Julian Day Observed Flows Modeled Flows 0 20000 40000 60000 80000 100000 120000 140000 160000 180000 200000 20 70 120 170 220 Vol um e ( cu bic me ter s) Julian Day Observed Flows Modeled Flows 54 Fig. 2.19 Validation of Mulberry Fork subwatershed 0 200 400 600 800 1000 1200 244 294 344 394 444 494 544 Flow (cm s) Julian Day Observed Flows Modeled Flows 0 50000 100000 150000 200000 250000 244 294 344 394 444 494 544 Flow (cm s) Julian Day Observed Modeled 55 Table 2.9 Calibration coefficients for the Mulberry Fork subwatershed SOIL LAYERS Layer 1 Layer 2 Layer 3 Prec Weight 0.95 Temp Lapse 4 Surface Roughness 0.2 Stream Roughness 0.115 Soil Thickness (cm) 15 18 20 Initial Moisture 0.5 0.5 0.5 Field Capacity 0.35 0.3 0.3 Saturation 0.45 0.4 0.35 Horizontal Conductivity (cm) 45000 45000 45000 Vertical Conductivity (cm) 22500 22500 22500 Root Distribution 0.75 0.1 0 Soil Density 1.3 1.3 1.3 Soil Tortuosity 10 10 10 2.5.6 Discussion of Calibration and Validation Results Calibration was successful for each of the four subwatersheds. The statistical analyses shown in Table 2.5 ranked the calibration for each subwatershed as a minimum of ?good? for both the calibration and validation periods, with the exception of the PBIAS for Blackwater Creek. The NSE coefficient was greater than 0.7 for both the calibration and validation periods for each subwatershed. The time periods for calibration and validation were chosen so that the model would cover both significant rain events in the spring, as well as the low-flow conditions through the late summer. During the calibration process, calibration coefficients were kept as consistent as was reasonable from one subwatershed to the other. The most sensitive coefficient is the precipitation weighting factor. The key coefficients of the soil parameters then controlled 56 how the system retained and discharged the moisture that reached the soil. The key parameters are the soil thickness, field capacity, saturation moisture, and conductivity values. These values are compared to literature values (Chen et al. 2001) in Table 2.10, using the data from the top layer of soil, which controls the majority of the watershed response. The final calibrated values for soil thickness match very well with the downloaded soil data given in 2.4.3 WARMF Simulation The simulation for WARMF was set to run on an hourly time step from January 2, 2010, to September 1, 2010. This time period was used as the calibration period to tweak the input parameters. After achieving the best calibration, the time period from September 2, 2010 to June 30, 2011 was used to validate the model calibration. This time frame allowed the model to be evaluated under a variety of conditions, including periods of heavy precipitation as well as the seasonal dry period in late summer. Table 2.3, given that the soils in each county cover a wide range of values. The thinnest soil is observed in Blount County (average of 16 cm), and the two subwatersheds from that area both achieved calibration using a top layer thickness of 15 cm. The calibration coefficients that are far out of the literature range are the conductivity values. The downloaded soil data suggest hydraulic conductivities of less than 1000 cm/day. However, calibrations were much worse when conductivities of that scale were attempted. The model needed conductivity values several of orders of magnitude larger. The model is very sensitive to input values of hydraulic conductivity, but relatively insensitive to the values for vertical conductivity. For calibration, this study used the standard practice of keeping the vertical conductivity at half the value of the horizontal conductivity (Dardashti 2010). 57 The high values for conductivity in the model are likely due to the very large catchment size caused by the 100-catchment limit in the program. The large catchments are assigned an average slope across the catchment. However, within each catchment, particularly within the rough, hilly terrain of this watershed, the slopes will vary considerably and the average slope will not be a good representation of the physical reality of the catchment. The overland flow portion of the flow routing in the model was likely insufficient, because the model gave very poor results when typical values were used for conductivity. The very high conductivity values compensated for the poor representation of the overland flow component. Table 2.10 Calibrated Parameters Compared to Literature Values Parameters Literature Range Calibration Range Precipitation Weighting Factor 0.5-1.4 0.8-0.95 Number of Soil Layers 1-5 3 Thickness of Soil Layers (cm) >0 15-22 Saturation Moisture 0.2-0.6 0.45 Field Capacity 0-0.4 0.35 Initial Moisture 0-0.6 0.4-0.5 Horizontal Conductivity (cm/day) >0 45000-190000 Vertical Conductivity (cm/day) >0 22500-95000 Even though very high values of hydraulic conductivity were input to overcome the poor modeling of overland flow, the calibration should not be discounted as being physically unrealistic. Exploring the model?s method of calculating hydraulic conductivity and the interaction between these coefficients demonstrates why these values give good results. WARMF uses Darcy?s Law (Equation 2.13) as the basis of its groundwater modeling (Chen et al. 2001). Qj is the lateral exfiltration rate, Khj is the horizontal hydraulic conductivity, S is 58 the slope of the subwatershed, W is the length of the catchment parallel to the receiving stream, and Zj is the thickness of the soil layer j. (2.13) The hydraulic conductivity (Khj) is a time-variable value in the model. At each time step, WARMF calculates the moisture of each soil layer using a mass balance approach. If the moisture of a layer is below the field capacity, then the hydraulic conductivity is zero, thus the exfiltration rate is zero. If the moisture of the layer is at saturation, then the hydraulic conductivity is taken to be the value input into the model, or K*. If the moisture in the soil layer is between the input field capacity and the saturation moisture, then the value of Khj to be used in Equation 2.12 is interpolated, based on Equation 2.14. In this equation, K*j is the input hydraulic conductivity value, ?j is the moisture of the soil, ?fc is the field capacity, and ?sj is the saturation capacity of the soil. (2.14) Thus the input values for hydraulic conductivity (K*), field capacity (?fc), and saturation moisture (?sj) are inextricably linked. If one value is adjusted, then the others are affected as well. For example, the absolute values for field capacity and saturation moisture are not as important as the difference between the values. If the two values are relatively close together, then exfiltration will either be nonexistent or the maximum possible. During calibration, it was observed that lowering field capacity and raising horizontal conductivity had similar effects, but were not equivalent. A lower value for field capacity increased the baseflow during low-flow periods, but the very high values for hydraulic conductivity were necessary to accurately model the falling limb for the event hydrographs. 59 2.6 Summary A hydrologic model for the Bankhead watershed was developed using WARMF. The model was developed using downloaded datasets through BASINS software as well as datasets created from other sources, including meteorological data from AWIS and streamflow data from the USGS. The model was calibrated using four USGS gage stations which divide the upper portion of the watershed into four subwatersheds: Blackwater Creek, Locust Fork, Lost Creek, and Mulberry Fork. The model was run on an hourly time step, and was calibrated from 01/02/2010 to 09/01/2010. The validation period was from 09/02/2010 to 06/29/2011. The calibration results were excellent, with NSE values above 0.75, and PBIAS values under ?3%, with the exception of Blackwater Creek (10%). As expected, validation results were not as exemplary, but were still satisfactory, with NSE values above 0.7. During the validation period, PBIAS errors were larger for Blackwater Creek (24.7%), Mulberry Fork (-13.1%), and Locust Fork (-15.8%). The graphical analyses agreed with the statistics, showing that the model is accurate to predict peak events if provided appropriate precipitation data. However, the PBIAS errors during the validation period suggest that the model could be a poor predictor if mass balance issues are important. In order to match observed flows, very high values for hydraulic conductivity were used. These high values were required in order to make up for the poor representation of the overland flow. The overland flow component of the flow routing was not adequate, because the large catchment sizes did not allow for realistic slope values. 60 Now that the watershed flow inputs have been modeled, the flows can be input into the EFDC hydrodynamic model, and the water balance utility can be used to develop the mass balance for the reservoir. 61 Chapter 3 Water Balance Analysis 3.1 Introduction/Background A hydrodynamic model was previously developed for Bankhead Reservoir, which is located on the Black Warrior River between Lewis Smith Dam and Bankhead Dam. This segment of the river is commonly called Bankhead Reservoir. The hydrodynamic model, developed using EFDC (Environmental Fluid Dynamics Code) (Hamrick 1992) had difficulties accurately modeling the hydrodynamics in the reservoir. The inaccuracy of the EFDC model is thought to be a result of a poorly quantified mass balance of the system. There are multiple small tributaries which flow into the Black Warrior River over the 75-mile domain which are not accounted for in the model. There is also known to be a leakage issue at Bankhead Dam. Therefore, a water balance approach was taken to attempt to improve the EFDC model. Certain inflows and outflows from the system are known, such as the inflow from Smith Dam, the discharge (non-leakage) through Bankhead Dam, the withdrawal intake for Plant Miller, and the drinking water intake for the Birmingham Water Works. Also, the four largest tributaries have flow data measured at USGS gage stations, but these stations are located significantly upstream from the reservoir, therefore the data do not represent the true flow entering the domain when measured flows from these USGS stations were used in EFDC model for Bankhead reservoir. 62 Other inflows are not easily quantified and require estimation or modeling. There is an unsteady, unknown leakage rate at Bankhead Dam as well as the water loss due to locking. Additionally, the primary unknown flow contribution is the watershed runoff, which includes the tributary inflows. In order to represent these flows, a watershed model was developed using WARMF (Watershed Analysis Risk Management Framework) (Chen et al. 2001). These inflows were then combined with the water balance tool contained in the CE- QUAL-W2 (Cole 2011) model to give a water balance that could be added into the EFDC model to represent the unknown flows. Proper quantification of the mass in the system will allow the model to represent the dynamics of the system more reasonably. 3.2 Water Balance Development Using WARMF 2010 and 2011 Low-flow Models The construction of a water balance for the reservoir hinged on an accurate hydrologic model to represent the flows contributed by the watershed. The previous WARMF model presented in Chapter 2, while it performed well over the calibration and validation time periods, was found to massively overpredict streamflows during the summer months when streamflow was at baseflow conditions during the dry summer months. The overprediction is shown in Fig. 3.4, demonstrating that the model overpredicted the total flows by as much as 400%. Because the purpose of this EFDC model was to model the reservoir through the summer months from Julian Day 235.5?307.5 (August 24, 2010?November 4, 2010), the previous ?high-flow? calibration of WARMF was not considered suitable to use for the water balance. Therefore, a low-flow hydrologic model targeted for these time periods was needed. 63 Fig. 3.1 2010 High-flow calibration of Blackwater Creek over low-flow period Fig. 3.2 2010 High-flow calibration of Locust Fork over low-flow period 0 1 2 3 4 5 6 7 8 9 10 235 245 255 265 275 285 295 305 Flo w (cm s) Julian Day Blackwater Creek Observed High-flow Calibration 0 10 20 30 40 50 60 70 80 90 100 235 245 255 265 275 285 295 305 Flo w (cm s) Julian Day Locust Fork Observed High-flow Calibration 64 Fig. 3.3 2010 High-flow calibration of Lost Creek over low-flow period Fig. 3.4 2010 High-flow calibration of Mulberry Fork over low-flow period 0 0.5 1 1.5 2 2.5 3 235 245 255 265 275 285 295 305 Flo w (cm s) Julian Day Lost Creek Observed High-flow Calibration 0 5 10 15 20 25 30 235 245 255 265 275 285 295 305 Flo w (cm s) Julian Day Mulberry Fork Observed High-flow Calibration 65 Fig. 3.5 Sum of the four USGS tributaries compared to the total USGS observed flows in 2010 3.2.1 Development of the WARMF Low-Flow Hydrologic Model The nature of the study allowed for the WARMF hydrologic model to be used in a non-standard way. Typically, a hydrologic model is developed, calibrated, and then validated over a separate time period. Theoretically, the validated model should then be able to be applied to other time frames with reasonable success. However, this traditional approach is not applicable for the low-flow study for two main reasons: 1) the time period of the study is small compared to typical hydrologic time scales 2) typical hydrologic models focus on peak flow evaluation at the expense of accuracy in base-flow situations. This study focuses on a 72-day time period. Hydrologic models are typically designed to simulate a time domain of years, rather than days. Small time frames negate the ?averaging? of errors over large time scales, which makes creating a statistically reasonable 0 20 40 60 80 100 230 240 250 260 270 280 290 300 310 Flo w (cm s) Julian Day WARMF High-Flow USGS Observed 66 model that can be applied to any short time scale an unrealistic task given the assumptions and lack of quality data common to hydrologic models. Hydrologic models are traditionally designed to model peak flows and frequently the same models poorly characterize base flows. In low-flow hydrology, the bulk of streamflow is gained from groundwater storage and can be highly dependent on the geology of the region. The understanding of low-flow generation mechanisms and which processes are important is still poorly understood (Smakhtin 2001). After a review of different types of hydrologic models, Staudinger et. al (2011) showed that the models that better predicted low- flow scenarios were not clearly influenced by one particular process within the model, but by the combination of different processes. Accurate hydrologic modeling for baseflow conditions is important, because water quality requirements are typically connected to low flows (Vogel and Fennessey 1995). WARMF?s physically-based, mechanistic design (Chen et al. 2005) should be suitable for yielding a physically meaningful and accurate low-flow calibration. The development of the WARMF low-flow model differed from a traditional approach in the method taken to achieve model calibration. The typical approach is to use split-sample time periods to separate data into calibration and validation time periods. The goal of the process was to model low-flow conditions for both 2010 and 2011, which yields itself to this split-sample method. However, 2010 and 2011 exhibited very different flow characteristics. Year 2010 exhibited a few small rain events, while 2011 was dominated by a single very large rain event. For example, in Locust Fork, the largest tributary of Bankhead Reservoir, the peak flow rate in 2010 in the time domain of interest was 17 cms. During the 67 same time period in 2011, the large rain event caused a peak flow of 720 cms. These drastically different flow conditions made development of a single low-flow model difficult. The final goal of the project was not to develop a single hydrologic model applicable to any condition, but to establish a reasonable model to represent the watershed flows over two distinct time periods. The nature of this goal allowed the development of two distinct low-flow hydrologic models, each calibrated for a particular year. Although this approach is not the standard approach of hydrologic modeling, for this project it is appropriate, with the acknowledgement that these models have little or no applicability to time periods outside of the range for which they were calibrated. To aid in developing reasonable low-flow hydrologic models, different evaluation criteria should be used than were previously employed. The Nash-Sutcliffe coefficient (NSE) and Percent Bias (PBIAS) were previously used in the calibration of the high-flow model. However, the Nash-Sutcliffe Coefficient uses squared errors, which magnify the error in peak flows at the expense of error in low flows. One proposed statistic to overcome this bias toward peak flows is the LNSE, which is essentially a modified Nash-Sutcliffe Efficiency that uses logged values of observed and predicted flows (Krause et al. 2005). Taking the natural logarithm of the flow values greatly reduces the statistical influence of the peak values, allowing low flows to have more influence. In addition to LNSE, PBIAS is still an appropriate measure of the models accuracy, because PBIAS is a direct measure of the volume of runoff, which is the primary concern of this project. Therefore, LNSE and PBIAS were the primary statistical measures used in conjunction with graphical analysis for the low- flow models. 68 3.2.2 WARMF Low-flow Model Setup Model output was needed over the EFDC time period of Julian Day 235.5-307.5, which is August 24 ? November 4 for both years. In order to give the models adequate spin- up time, the model was run from Julian Day 151 (June 1). 3.2.2.1 Recalibration of WARMF Model for Low-flow Conditions Calibrating the WARMF model for low-flow conditions involved re-examination of which parameters control the watershed response during baseflow conditions. The hydraulic conductivity, which controls Darcy?s Law (Equation 2.13), is the principle parameter for adjusting the baseflow modeling. Within the hydraulic conductivity are the parameters of field capacity and saturation moisture (Equation 2.14). Therefore, the parameters shown in Table 3.1 were adjusted to match the baseflow conditions for the 2010 and 2011 low- flow WARMF models. Table 3.1 Calibration Parameters for WARMF Low-flow Calibration Blackwater Locust Lost Mulberry Parameter 2010 2011 2010 2011 2010 2011 2010 2011 Precipitation Weight 0.8 1 0.82 1 0.95 1 0.95 1 Temp Lapse 3.5 0 0 0 5 0 4 0 Soil Thickness (cm) ? Layer 1 20 22 20 22 20 20 20 22 Soil Thickness (cm) ? Layer 2 22 25 22 25 35 22 25 Field Capacity 0.15 0.3 0.17 0.25 0.17 0.3 0.35 0.3 Saturation Moisture 0.45 0.45 0.45 0.45 0.3 0.45 0.45 0.45 Horizontal Conductivity (cm/day) 1300 14000 1300 10000 1500 14000 800 8000 3.2.2.2 Analysis of 2010 and 2011 Low-flow Calibrations The statistical results for both the 2010 and 2011 low-flow models are shown in 69 Table 3.2. Detailed discussion of the results is in the following sections. In the 2010 model, only Lost Creek meets the established NSE criteria of 0.5. However, when the LNSE is considered in order to better evaluate the model?s performance modeling baseflow, the results are comparatively improved. LNSE has not be sufficiently used in literature for an acceptable minimum value to be established (Krause et al. 2005), but it is assumed that overall, the evaluation criteria for NSE and LNSE would be similar. Three of the four sites exhibit very good bias control under ?10%, with Mulberry Fork as an outlier. However, an explanation for this variation will be presented. The 2011 model showed better performance when evaluated using the NSE and LNSE parameters. The NSE was remarkably high for this model, due to the nature of the watershed flows for this model. Unlike the consistent low flows in 2010, the 2011 time period is dominated by a single large rain event. Thus the NSE will become dominated by whether the models predict this large peak accurately. The high NSE values show that the model does capture this peak well, and higher LNSE values show that the lower flows are also fairly well-represented with each site having an NSE > 0.5. The PBIAS values are very good, with each value under 10%. Table 3.2 Statistical results for the 2010 and 2011 Low-flow models 2010 2011 NSE PBIAS (%) LNSE NSE PBIAS (%) LNSE Blackwater Creek 0.16 -5.24 0.48 0.84 0.31 0.63 Locust Fork 0.09 -7.22 0.37 0.89 7.25 0.55 Lost Creek 0.86 3.19 0.78 0.89 6.57 0.70 Mulberry Fork 0.45 -19.64 0.81 0.91 2.37 0.77 BLACKWATER CREEK 70 In 2010, the Blackwater Creek hydrograph (Fig. 3.6) shows that the model misses both the magnitude and timing of the small peak events. These events have a maximum magnitude of 3.0 m3/sec, which is a small amount for the size of the watershed. These small discrepancies are likely due to input data. Prediction for small changes in flow requires very fine input data. The base flow conditions are matched fairly well, which is evident in the cumulative curve and the PBIAS value of -5.23%. In 2011 (Fig. 3.7), the model matches the peak flow fairly well. The cumulative curve shows that the model tracks with the observed volume well, which indicates an accurate prediction of baseflow conditions. LOCUST FORK For the 2010 calibration, the Locust Fork (Fig. 3.8) models the small peaks more effectively than the Blackwater Creek site, but predicts two larger peaks where the observed data do not have a peak. This difference is likely due to rainfall input data, because the other peaks are matched more accurately. The poor temporal match of the peaks combined with the two false predicted peaks results in a low NSE value. Inspection of the cumulative discharge curve shows that with the exception of the error caused by the false peak at Julian Day 270, the cumulative discharge is predicted almost perfectly, because the baseflow is represented very well. The 2011 model (Fig. 3.9) predicts the peak flows fairly well with the exception of underestimating the falling limb of the hydrographs. The cumulative discharge curve shows that the underestimation is largely contained in these two events and that the model represents the base flow fairly well. 71 LOST CREEK In 2010, Lost Creek is essentially a baseflow-only model, with the only peak a small (~1 cms) increase in flow (Fig. 3.10). Also of note, up to Julian Day 240, the observed data appear to be flawed, because the observed hydrograph appears to be a step function. The cumulative distribution confirms that the model accurately describes the baseflow in the subwatershed. In 2011 (Fig. 3.11), the trends observed in Blackwater Creek and Locust Fork watersheds continued. Baseflow is represented well, but the total volume contributed by the large rain event was underestimated by poorly representing the falling limb of the peak hydrograph. MULBERRY FORK For the 2010 model (Fig. 3.12), the cumulative curve shows that the volume of the largest rain event is overpredicted. However, after that overprediction, the predicted cumulative curve tracks the observed volume very well, which indicates that the baseflow is modeled very well. The high PBIAS of -19.6% is caused solely by this single overestimation. In the 2011 model (Fig. 3.13), the peak flow was modeled almost perfectly, reflected by the high NSE value of 0.91, which is confirmed by the hydrograph and the cumulative discharge curve. At Julian Day 264, there is a small observed peak which is not reflected in the model, which causes the difference in the observed and modeled cumulative volume curves after that point. 72 Fig. 3.6 Calibration of Blackwater Creek for 2010 low-flow model 73 Fig. 3.7 Calibration of Blackwater Creek for 2011 low-flow model 74 Fig. 3.8 Calibration of Locust Fork for 2010 low-flow model 75 Fig. 3.9 Calibration of Locust Fork for 2011 low-flow model 76 Fig. 3.10 Calibration of Lost Creek for 2010 low-flow model 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 2 . 5 3 . 0 3 . 5 4 . 0 4 . 5 5 . 0 230 240 250 260 270 280 290 300 310 D i s c ha r g e ( c m s ) J u l i a n Da y Di s c h a r ge o ve r E F DC T i m e Do m a i n ( 2 0 1 0 ) W A R M F L o s t U SG S L o s t 0 100 200 300 400 500 600 700 800 230 240 250 260 270 280 290 300 310 D i s c ha r g e (c ub i c m e t e r s ) J ul i a n D a y C u m u l a t i ve Di s ch a r ge o ve r E F DC T i m e Do m a i n ( 2 0 1 0 ) W A R M F C u m u la t iv e L o s t U SG S C u m u la t iv e L o s t 77 Fig. 3.11 Calibration of Lost Creek for 2011 low-flow model 78 Fig. 3.12 Calibration of Mulberry Fork for 2010 low-flow model 79 Fig. 3.13 Calibration of Mulberry Fork for 2011 low-flow model 80 3.3 Background of CE-QUAL and EFDC Hydrodynamic Models After development of reasonable watershed input through the WARMF low-flow models, the water balance of the reservoir could be performed with the CE-QUAL-W2 water a balance utility and finally the hydrodynamic processes within the reservoir could be simulated using EFDC. Hydrodynamic and water quality models have proliferated in recent years due to availability of computational power and developments in numerical modeling (Batick 2011). The models range from simplistic 1-D models to sophisticated 3-D models. EFDC (Environmental Fluid Dynamics Code) (Hamrick 1992) and CE-QUAL-W2 (Cole 2011) are both hydrodynamic and water quality models that rely on the conservation of mass and momentum to model hydrodynamics, temperature, and other constituents. This study uses a pre-existing EFDC 3-D model, so a separate full CE-QUAL-W2 2-D model was not used. Only the water balance utility which is packaged with CE-QUAL-W2 was used. 3.3.1 CE-QUAL-W2 CE-QUAL-W2 is a 2-dimensional, laterally averaged hydrodynamic water quality model (Cole and Buchak 1995). It is continuously updated and distributed through Portland State University. The most recent update, Version 3.7, includes support for comma-delimited input file formats and a postprocessor (Cole 2011). CE-QUAL-W2 is designed to be able to model any number of reservoirs, branches, and tributaries simultaneously. CE-QUAL-W2 predicts water surface elevations, velocities, temperatures, and many constituent concentrations. It has been successfully applied to hundreds of reservoirs around the world, 81 and is highly adaptable to most natural conditions as long as the model?s limitations are considered. The model does not consider Coriolis force, so very large water bodies may not be appropriate to be modeled. Additionally, the lateral-averaging scheme limits the model to water bodies where variation of simulated parameters in the third dimension is negligible, which generally means the model is appropriate for dendritic systems (relative narrow reservoirs in comparison to the reservoir length). The model also assumes that the length of the reservoir is much greater than the depth of the reservoir, which allows simplification of the z-momentum equation (Batick 2011). The model features an auto-stepping algorithm that dynamically determines the maximum time step which yields stable results. The model includes a basic user interface that allows the user to set and change model parameters. However, there is no interface to aid in developing the bathymetry file, other input text files, or to view results. 3.3.2 CE-QUAL-W2 Water Balance Utility Included with the model is a water balance utility designed to be the first step in calibrating the model. First, the model is run without a water balance to give an initial water surface elevation, typically taken at the most downstream point. The utility then compares the modeled surface elevations to the observed elevations and outputs a flow file that theoretically should bring the reservoir into balance. The water balance process typically requires several iterations, with the resulting flow file being changed until the modeled elevations match the observed elevations. Thus the required files for the water balance utility are as follows: 82 W2 control file (w2con.npt) ? Contains the start and end date for balance analysis W2 bathymetry file (bth.npt) ? Contains volume-elevation information for calculating amount of volume needed to achieve balance Time series of modeled elevation (tser_1_seg.opt) ? Modeled water surface elevation at comparison point (typically most downstream point) Observed water surface elevations (el_obs.npt) ? Observed elevations to compare modeled elevations Flow balance output file (qwb.opt) ? Required flow to add/remove into the water body to achieve mass balance Although only the water balance tool of CE-QUAL-W2 and not the full hydrodynamic model was to be used, the development of the reservoir bathymetry was still required. HEC-RAS and Corps of Engineers bathymetric data were available for the reservoir, so these data were converted into a CE-QUAL-W2 bathymetry file (Cole 2011). To ensure that the CE-QUAL-W2 and EFDC bathymetries were equivalent, and thus the water balance approach could be valid, the volume-elevation curves were compared (Fig. 3.14). The CE-QUAL model at Bankhead and Plant Gorgas is shown in Fig. 3.15 and Fig. 3.16. 83 Fig. 3.14 Comparison of EFDC and CE-QUAL bathymetries for the whole reservoir Fig. 3.15 View of CE-QUAL model near Bankhead Dam 0 50 100 150 200 250 300 40 45 50 55 60 65 70 75 80 Vo lum e ( 10 6 c ub ic m eter s) Elevation (m) Volume Comparison EFDC CE-QUAL 84 Fig. 3.16 View of CE-QUAL model near Plant Gorgas 3.3.3 EFDC (Environmental Fluid Dynamics Code) Hydrodynamic Model EFDC is a 3-D hydrodynamic model that simulates temperature, water quality, sediment transport, and contaminant transport by solving 3-D equations for continuity, momentum, and the free water surface (Hamrick 1992). EFDC has been applied to numerous study areas, including the James River and New York estuaries (Hamrick et al. 1995), Lake Okeechobee (Jin et al. 2002), thermal discharge into Conowingo Pond in New York (Hamrick and Mills 2000), Lake Tenkiller in Oklahoma (Ji et al. 2003), and Morro Bay, California (Ji et al. 2001) among others. EFDC is a public-domain code that was originally developed at the Virginia Institute of Marine Science. The U.S. EPA contracted Tetra Tech, Inc. to develop EFDC-Hydro, 85 which is the publicly available model. This study uses EFDC-Explorer, which is a GUI- equipped modified version of EFDC from Dynamic Solutions-International, LLC. 3.4 Development of EFDC Model to Incorporate Water Balance After development of the 2010 and 2011 WARMF hydrologic models, the watershed flows were then exported for inclusion in a total water balance model of Bankhead Reservoir. A total of 8 tributary inflows (including the four gaged streams used for calibration) and 11 catchment inflows were represented and are shown in Fig. 3.18. The streamflow gages on the four major tributaries are located significantly upstream from the reservoir, so rather than using the gage values in the reservoir model, the predicted WARMF output of the streams at the reservoir was used. Fig. 3.17 shows the difference between the WARMF modeled flows at the USGS gage and at the intersection with Bankhead Reservoir for the 2010 year of Locust Fork. The inflows from the 11 catchments are the predicted surface runoff from the catchments located directly on the reservoir, whose runoff is not included in another tributary. The horizontal grid for the EFDC model contained 6,974 horizontal cells, which each had 10 layers for a total of 69,740 computational cells. The DX values ranged from 9.5 m to 189.8 m with an average DX of 25.7 m. The DY values ranged from 10.0 m to 277.1 m, with an average DY of 100.9 m. 86 Fig. 3.17 Comparison of 2010 WARMF modeled flows of Locust Fork at the USGS gage and at the intersection with Bankhead Reservoir 0 5 10 15 20 25 30 235 245 255 265 275 285 295 305 Flo w (cm s) Julian Day Flow at USGS gage Flow at Bankhead Reservoir 87 Fig. 3.18 Locations of boundary condition tributaries in the EFDC model domain. Major locations of interest are labeled. Other points are minor tributaries or the pour point of a subcatchment located on the reservoir. 3.4.1 EFDC Boundary Condition Setup The EFDC model was modified to include the tributary inflows as boundary conditions. After modification to the EFDC model, the model can be summarized as follows: Upstream Boundary Condition ? Smith Dam inflow Downstream Boundary Condition ? Bankhead Dam outflow Water Balance Boundary Condition ? A set of three cells near Bankhead Dam were assigned to manage the future water balance flows 88 Withdrawal/Return Condition ? Set to model the effect of Plant Gorgas on the reservoir Previous Boundary Conditions ? Sipsey drinking water withdrawal, Plant Miller Intake withdrawal, Mulberry Fork drinking water intake Watershed Input Boundary conditions ? The previously described 8 tributaries and 11 catchment inflows were assigned to the reservoir at their point of confluence. For catchments, the flows were added at the downstream point of the catchment. Under a traditional usage of the CE-QUAL-W2 water balance utility, the outputted flow balance file (qwb.opt) would be added back into the reservoir as a distributed tributary, which distributes the flow equally into all segments of the water body. This design is to simulate groundwater supply/removal, evaporation, or some other consistent mechanism that inputs or removes flow from the length of the water body. For this application, the water balance was intended to rectify inputs from the watershed, but was also intended to roughly model the effect of a known leakage issue at Bankhead Dam and the loss of water through locking. Because EFDC does not have a distributed tributary function, the entire mass balance flow condition was applied near the downstream boundary condition of Bankhead Dam. The downstream boundary conditions as set up in EFDC are shown in Fig. 3.20. The time-series flow files of Smith Dam and Bankhead Dam for 2010 and 2011 are shown from Fig. 3.21 to Fig. 3.24. EFDC was also set up with the required atmospheric forcing data, including wind speed and direction, air temperature, atmospheric pressure, relative humidity, precipitation, cloud cover, and solar radiation. The hourly atmospheric data were obtained from the AWIS 89 site in Cleveland, AL, which was the closest site that could provide solar radiation data. Atmospheric data for 2010 and 2011 are shown from Fig. 3.25 to Fig. 3.34. Fig. 3.19 Upstream boundary of EFDC model Fig. 3.20 Downstream boundary conditions of EFDC model 90 Fig. 3.21 2010 upstream flow boundary at Smith Dam Fig. 3.22 2011 upstream flow boundary at Smith Dam 91 Fig. 3.23 2010 downstream flow boundary at Bankhead Dam Fig. 3.24 2011 downstream flow boundary at Bankhead Dam 92 Fig. 3.25 2010 air temperature from Birmingham Airport Fig. 3.26 2011 air temperature from Birmingham Airport 0 5 10 15 20 25 30 35 40 45 50 235 245 255 265 275 285 295 305 De gre es Ce lsius Julian Day 93 Fig. 3.27 2010 relative humidity from Birmingham Airport Fig. 3.28 2011 relative humidity from Birmingham Airport 94 Fig. 3.29 2010 rainfall data from Birmingham Airport Fig. 3.30 2011 rainfall data from Birmingham Airport 95 Fig. 3.31 2010 solar radiation data from Birmingham Airport Fig. 3.32 2011 solar radiation data from Birmingham Airport 96 Fig. 3.33 2010 wind speed data from Birmingham Airport Fig. 3.34 2011 wind speed data from Birmingham Airport 97 Plant Gorgas was modeled through withdrawal/return boundary conditions (Fig. 3.35). Water is drawn into the plant?s intake canal through an underwater weir on the river. From the intake canal, there are two cooling water intake structures (CWIS), one which provides water for units 6 and 7, and one with provides water for units 8, 9, and 10. In reality, the intake for units 8, 9, and 10 is routed through a pipe underneath the discharge canal (Baker?s Creek) into a small pond, but for simplicity, the model was set up to pull intake water directly from the intake canal. The aerial view showing the intake/discharge areas is given in Fig. 3.36. Fig. 3.35 Withdrawal/return boundary conditions simulating Plant Gorgas Fig. 3.36 Intake and discharge paths for Plant Gorgas 98 3.5 EFDC Model Methodology Separate EFDC models for 2010 and 2011 were run from Julian Day 235.5-307.5 (August 24?November 4). The initial water surface elevations were set to the observed downstream elevations at the start time, which were 77.605 m in 2010 and 77.547 m in 2011. Initial temperatures for the reservoir were taken from observed temperature profiles. Using observed profiles as initial temperature conditions allows the model to have a much shorter spin-up time. If a constant temperature was used throughout the reservoir, then deep, stratified layers may not be modeled correctly, or an unnecessarily large spin-up time would be required. In order to evaluate the effect of the water balance approach, two models were run for each year. A ?default? model that did not include any water balance or WARMF input, but did include the measured USGS data, was ran initially. These models are called the ?2010 NO BALANCE? and ?2011 NO BALANCE? models. The water surface elevations from the downstream boundary at Bankhead Dam were then exported and formatted to function inside the CE-QUAL water balance tool. The water balance tool analyzed the modeled data, compared it to the observed elevation data, and output a flow file that theoretically would bring the reservoir into balance. This flow file was output on an hourly timestep and input into EFDC as a boundary condition near Bankhead Dam as previously shown. The output flow balance file can be output at various time intervals and can also be averaged using a running average over any time period. After comparing various combinations of skip intervals and averaging intervals, the hourly time interval with an hourly running average was used. Using longer periods of running average yields a smoother 99 and more consistent mass balance curve; however, this method does not match water surface elevations as consistently as a mass balance with an hourly timestep. After re-running the models with the final flow balance, the modeled downstream elevations were brought into better agreement but still had significant deviations. The modeled elevations were again exported to the water balance utility and an updated flow file created. This iterative process was followed until good agreement was reached between the observed elevations and the modeled elevations. Several iterations were required, partially because neither initial model completed the whole time domain before crashing. The 2011 model reached a balance after five iterations. The 2010 model required more iterations because different time steps and averaging periods were tested during the process. For both models, the average flow rate followed the observational loss of flow through Bankhead Dam. The final balance used in the EFDC model was computed on an hourly time scale which resulted in fairly large, rapid fluctuations from negative to positive flows to keep the reservoir in balance on an hour-to-hour basis. The time-series mass balance flow for the 2010 BALANCED model is shown in Fig. 3.37. and the 2011 balance is shown in Fig. 3.38. 100 Fig. 3.37 2010 Mass balance flow boundary added near Bankhead Dam Fig. 3.38 2011 Mass balance flow boundary near Bankhead Dam 101 3.6 EFDC Model Results and Analysis The final downstream water surface elevation at Bankhead Dam agreed well with the observed elevations. This result was naturally expected, because the methodology was designed to optimize this relationship. In order to evaluate the model?s performance throughout the rest of the water body, the water surface elevation at the USGS station at Cordova was used. Fig. 3.39 Location of Cordova USGS gage 3.6.1 2010 BALANCED EFDC Model Fig. 3.40 shows the water surface elevation at Bankhead Dam compares the 2010 NO BALANCE model, the 2010 BALANCED model, and the observed elevations. The 2010 NO BALANCE model is predictably poor, with the elevations continuously rising, while the observed elevations stay relatively constant. This continuous rise clearly demonstrates that there is flow leaving the real system that is not correctly represented in the model, which reinforces the theory that significant water is lost through leakage and locking from Bankhead Dam. The water surface elevation rise caused the EFDC model to crash in Julian Day 298.Fig. 3.41 is a closer examination of the difference between the final balanced model and the observed elevations at Bankhead. The modeled elevations track the observed elevations very well, which should be the case, considering that the reservoir was balanced 102 based on this datapoint. The average absolute error in the modeled data is only 2.3 cm. A statistical summary showing average and average absolute errors for the two 2010 runs is shown in 103 Table 3.3. Fig. 3.42 shows the water surface elevations at the Cordova USGS station for both the 2010 NO BALANCE and 2010 BALANCED runs. From Julian Day 286-296, the gage appears to have malfunctioned. There is no physical explanation for this rise in surface elevation, and the sudden drop at day 286 suggests that this is the time when the gage resumed functioning properly. For the data analysis in 104 Table 3.3, this time period was removed from the analysis. In Fig. 3.42, the Cordova elevation from the 2010 NO BALANCE model again rises continuously, as would be expected from the previous analysis. Fig. 3.43 more clearly demonstrates the 2010 BALANCED model?s behavior at the Cordova gage. The model follows the observed elevations very well, with an average absolute error of only 5.3 cm. 105 Table 3.3 Statistical summary of elevation differences between the 2010 NO BALANCE and 2010 BALANCED models Cordova Elevation (m) Bankhead Elevation (m) BALANCED NO BALANCE BALANCED NO BALANCE 2010 Avg Error -0.031 0.951 0.018 1.365 2010 Abs Error 0.053 0.955 0.023 1.366 2011 Avg Error 0.029 -0.174 -0.005 -2.576 2011 Abs Error 0.150 0.286 0.021 2.710 3.6.2 2011 BALANCED EFDC Model The 2011 calibration was performed following the same methodology as the 2010 calibration. In Fig. 3.44, the 2010 NO BALANCE model did not perform well, and crashed at Julian day 250. Examining the time series of modeled water surface elevations shows that the model crash is caused by the large flows due to the storm event. At Bankhead, the large volume of water being removed through the spillways is not being replaced by watershed input as would be the case in reality. The observed elevations remain fairly stable, but the modeled elevations have a massive drop of 15 m. After the balance procedure, the elevations matched well, as shown in Fig. 3.45, with an average absolute error of 2.1 cm. At Cordova, shown in Fig. 3.46, the water surface elevation spikes during the flood as expected. The 2011 NO BALANCE model does not increase as high as the observed elevations. Even though the 2011 NO BALANCE model includes the large releases from Smith Dam, the lack of representing all of the tributary inflows causes an underestimation in the model. The 2011 BALANCED model matches the observed elevations very well, with an average absolute error of 15 cm. 106 3.6.3 Summary of EFDC Modeled Elevations For both the 2010 and 2011 models, the water balance flow file yielded a model that very accurately represented the water surface elevation throughout the reservoir. This result was expected because of the nature of the optimization procedure. However, the next step of the analysis is to determine whether or not this flow balance file represents what is truly happening in the reservoir, or whether it is simply forcing a type of pressure boundary condition. The negative mass flows are physically explained by leakage from the dam and loss through locking. The positive mass flow rates are not as easily explained. In general, the positive values are a result of the small hourly timestep of the water balance procedure. The procedure repeatedly over-corrects itself, which results in the mass flow swinging back and forth across a mean negative flow. As previously discussed, a longer timestep can be used in the flow balance, but it does not as accurately represent the water surface elevation on a finely resolved timescale of minute data. Some mass flows can be physically explained for several reasons. The ash pond discharge of Plant Gorgas contributes a steady small flow to the reservoir. The WARMF low-flow models tend to underestimate the tributary flow. Additionally, there may be groundwater recharge along the reservoir that is unaccounted for by the surface-layer WARMF model. In order to specifically evaluate whether the flow boundary condition more accurately represents the hydrodynamics of the reservoir, then other constituents would need to be evaluated, along with velocities at cross sections of the waterbody. Due to a lack of data, a 107 conceptual discussion is the limit of what can be currently undertaken. However, future study possibilities will be discussed in the next section. 108 Fig. 3.40 2010 comparison of modeled and observed water surface elevations at Bankhead Dam 109 Fig. 3.41 2010 comparison of final balanced water surface elevation to the observed elevation at Bankhead Dam 110 Fig. 3.42 2010 comparison of modeled and observed water surface elevations at Cordova USGS station 111 Fig. 3.43 2010 comparison of final balanced water surface elevation to the observed elevations at Cordova USGS station 112 Fig. 3.44 2011 comparison of modeled and observed water surface elevations at Bankhead Dam 113 Fig. 3.45 2011 comparison of final balanced water surface elevation to the observed elevation at Bankhead Dam 114 Fig. 3.46 2011 comparison of modeled and observed water surface elevations at Cordova USGS station 115 Fig. 3.47 2011 comparison of final balanced water surface elevation to the observed elevation at Cordova USGS station 116 Chapter 4 Summary and Future Study 4.1 Summary of Study A hydrologic watershed model was developed for constructing a water balance on Bankhead Reservoir and providing watershed (stream) inflow data for a hydrodynamic model. Even though water levels in the reservoir are held fairly constant, the processes within the reservoir are highly dynamic. The timing and amount (discharge) of the upstream and downstream controlled releases create seiches and density currents. The presence of a power plant in the middle of the reservoir with once-through cooling complicates the dynamics of the system. The hydrodynamic model EFDC was used to represent the system. However, the model encountered difficulties that were believed to be a result of improper mass balance. The tributary inflows can contribute a significant volume to the reservoir, and the downstream control, Bankhead Dam, is known to have water loss through leakage and locking. The initial WARMF model was developed by the typical hydrologic modeling processes, and calibrated and validated using a split-sample time period. Although the initial calibration was good for the traditional application to high flow periods, the calibration significantly overpredicted baseflow and low-flow conditions. Because the hydrodynamic model was designed for use during the low-flow warm summer months, the initial high-flow calibration of WARMF was not appropriate to use as the watershed model for the mass balance study. 117 The model was to be run in the summer months (August-November) of 2010 and 2011. Therefore optimized WARMF models were developed for each summer period. Although this approach is not standard for hydrologic modeling, the approach was valid because the models were not assumed to be valid for any other time period than the time period for which they were optimized. Four USGS gage stations within the watershed provided the data by which the models were optimized to predict flows across the watershed. The time periods of 2010 and 2011 were very different hydrologic conditions. 2010 was a very dry year, with peak observed flows at the large tributary of Locust Fork only reaching 17 cms. 2011 was characterized by a single very large event, which caused peak flows of 720 cms in Locust Fork. The 2010 model was an acceptable calibration, with log NSE values above 0.37, and the 2011 model calibrated better, with all four log NSE values above 0.55. Most importantly, the flow biases were small, with the discrepancies being explained by input data. The water balance was calculated using the water balance utility in the CE-QUAL- W2 hydrodynamic model. The full 2-D CE-QUAL model was not used in favor of the 3-D EFDC model. The output of the water balance utility was incorporated into EFDC near the downstream boundary condition, and the lake water balance was determined based upon observed water surface elevations at Bankhead Dam. As expected, after iteratively running the water balance process, the modeled surface elevations in EFDC closely matched the observed elevations at Bankhead Dam. Bringing the reservoir into the proper mass balance should then allow for more complicated hydrodynamic and water quality processes to be simulated within the reservoir in the future. 118 4.2 Recommendations for Future Study The mass balance could be evaluated to determine if a more accurate method of determining watershed input could be used. Hydrologic models are neither accurate nor precise predictors of tributary flows. A range of options could be considered, from a simple drainage area ratio to complex hydrologic models coupled to full groundwater models. More gages and monitoring stations would be required to increase the accuracy of the models. The weaknesses observed in the developed WARMF models could be addressed through the application of a different model, or use of the full version of WARMF that does not have a limit on catchment size. The large catchment size used in this study caused issues by lumping diverse land conditions into catchments. The issues encountered in the model with the poor representation of overland flow and the excessively high values for hydraulic conductivity would probably be avoided with smaller catchment sizes. Now that the reservoir has been brought into a mass balance condition, there is potential for extensive future study. The reservoir has interesting hydrodynamic situations. The cold water release from Smith Dam causes a plunging inflow into the reservoir. The deep-layer intake for Plant Gorgas brings in the deep cold water, and surface-level drinking water intakes upstream of the plant cause a split-layer flow dynamic, where the cold deep water flows downstream, and the warmer surface water flows upstream. Any number of constituents (temperature, DO, sediment) could be studied to evaluate their impact on the reservoir under a variety of management scenarios. Creating a physically meaningful water balance is the key to unlocking the potential the reservoir for future study. 119 References Ambrose, R., Wool, T., and Martin, J. (1993). "WASP 5, The Water Quality Analysis Simulation Program Version 5.00." ASCI Corporation, Athens, Georgia. Arnold, J. G., and Soil, G. (1994). SWAT (Soil and Water Assessment Tool), Grassland, Soil and Water Research Laboratory, USDA, Agricultural Research Service. Batick, B. M. (2011). "Modelign Temperature and Dissolved Oxygen in the Cheatham Reservoir with CE-QUAL-W2." Vanderbilt University. Beasley, D. B., Huggins, L. F., and Monke, E. J. (1980). "ANSWERS: A model for watershed planning." Trans. of the ASAE, 23(4), 938-944. Bingner, R. L. (1996). "RUNOFF SIMULATED FROM GOODWIN CREEK WATERSHED USING SWAT." Transactions of the ASAE, 39(1), 85-90. Buck, A. L. (1981). "New Equations for Computing Vapor Pressure and Enhancement Factor." Journal of APplied Meteorology, 20, 1527-1532. Chapra, S. C. (1997). Surface water-quality modeling, McGraw-Hill New York. Chen, C. H., J.; Tsai, W. (2006). "Enhancement of Watershed Analysis Risk Management Framework (WARMF) for Mercury Watershed Management and Total Maximum Daily Loads (TMDLs).", EPRI, Palot Alto, CA and Minnesota Power, Duluth, MN. Chen, C. W., Gherini, S. A., Hudson, R., and Dean, J. (1983). "The Integrated Lake- Watershed Acidification Study. Volume 1: Model principles and application 120 procedures." Electric Power Research Institute report EA-3221, Palo Alto, California, USA. Chen, C. W., Herr, J., Weintraub, L. H. Z., Goldstein, R. A., Herd, R., and Brown, J. M. "Framework to calculate TMDL of acid mine drainage for Cheat River Basin in West Virginia." Proc., Watershed Management and Operations Management 2000, June 20, 2000 - June 24, 2000, American Society of Civil Engineers. Chen, C. W., and Herr, J. W. (2002). "Comparison of BASINS and WARMF models: Mica Creek watershed." Electric Power Research Institute, Palo Alto, CA. Chen, C. W., Herr, J. W., Goldstein, R. A., Ice, G., and Cundy, T. (2005). "Retrospective comparison of watershed analysis risk management framework and hydrologic simulation program Fortran applications to Mica Creek watershed." Journal of Environmental Engineering, 131, 1277. Chen, C. W., Herr, J. W., and Weintraub, L. (2001). "Watershed Analysis Risk Management Framework (WARMF): Update One?A decision support system for watershed analysis and total maximum daily load calculation, allocation and implementation. Publication No. 1005181." Electric Power Research Institute, Palo Alto, California. Chen, C. W., Loeb, C., and Herr, J. W. (2000). "Adaptation of WARMF to Calculate TMDL for Chartiers Creek Watershed in Pennsylvania." Final Rep. to USEPA Region 3, Philadelphia, PA and The Chartiers Creek Watershed TMDL Stakeholder Group, Pittsburgh, PA. Chu, T. W., and Shimohammadi, A. (2004). "EVALUATION OF THE SWAT MODEL?S HYDROLOGY COMPONENT IN THE PIEDMONT PHYSIOGRAPHIC REGION OF MARYLAND." Transactions of the ASAE, 47(4), 1057-1073. 121 Cole, R. W., and Buchak, E. M. (1995). "CE-QUAL-W2: A two dimensional, laterally averaged, hydrodynamic and water quality model. Version 2.0." U.S. Army Engineer Waterways Experiment Station, Vicksburg, Mississippi. Cole, T. M. W., Scott A. (2011). "CE-QUAL-W2: A Two-Dimensional, Laterally Averaged, Hydrodynamic and Water Quality Model, Version 3.71 User Manual."US Army Engineer Waterways Experiment Station, Vicksburg, MS. Dardashti, S. (2010). "Modeling Hydrology and Nitrogen Fate and Transport in a Tile- Drained Agricultural Watershed in a Cold Region." McGill University. Dayyani, S., Prasher, S. O., Madani, A., and Madramootoo, C. A. (2012). "Impact of climate change on the hydrology and nitrogen pollution in a tile-drained agricultural watershed in Eastern Canada." Transactions of the ASABE, 55(2), 389-401. Donigian, A. S., Jr., Imhoff, J. C., Bicknell, B. R., and Kittle, J. L., Jr. (1984). "Application guide for Hydrological Simulation Program--Fortran (HSPF)." EPA-600/3-84-065, U.S. Environmental Protection Agency, Environmental Research Laboratory, Athens, GA. Gao, C.-c., Zhu, M.-f., and Qin, H.-x. (2008). "Adaptability of Hargraeves Formula in Estimating Reference Evapotranspiration in Xinjiang Region [J]." Journal of Irrigation and Drainage, 3, 017. Geza, M., and McCray, J. E. "Modeling the effect of population growth on stream nutrient concentration in turkey creek watershed using the WARMF model." Proc., 11th National Symposium on Individual and Small Community Sewage Systems, October 20, 2007 - October 24, 2007, American Society of Agricultural and Biological Engineers. 122 Geza, M., and McCray, J. E. (2008). "Effects of soil data resolution on SWAT model stream flow and water quality predictions." J. Environ. Manag., 88(3), 393-406. Geza, M. M., John E (2009). "A User Guide for Modeling Watershed-Scale Impacts of Onsite Wastewater Systems: Case Studies of Impacts of Onsite Systems in Turkey Creek Watershed, Colorado." Colorado School of Mines, Environmental Science and Engineering Division. Gherini, S., Mok, L., Hudson, R., Davis, G., Chen, C., and Goldstein, R. (1985). "The ILWAS model: formulation and application." Water, Air, & Soil Pollution, 26(4), 425-459. Gollamudi, A. (2006). "Hydrological and Water Quality Modeling of Agricultural Fields in Quebec." M.Sc., McGill University. Gupta, H., Sorooshian, S., and Yapo, P. (1999). "Status of Automatic Calibration for Hydrologic Models: Comparison with Multilevel Expert Calibration." Journal of Hydrologic Engineering, 4(2), 135-143. Hamrick, J., Kuo, A., and Shen, J. (1995). "Mixing and dilution of the Surrey Nuclear Power Plant cooling water discharge into the James River." A report to Virginia Power Company, The College of William and Mary, Virginia Institute of Marine Science, Gloucester Point, VA. Hamrick, J. M. (1992). "A three-dimensional environmental fluid dynamics computer code: Theoretical and computational aspects." Speical Report 317, The College of William and Mary, Virginia Institute of Marine Science, Gloucester Point, VA, 63. 123 Hamrick, J. M., and Mills, W. B. (2000). "Analysis of water temperatures in Conowingo Pond as influenced by the Peach Bottom atomic power plant thermal discharge." Environmental Science & Policy, 3(Supplement 1), 197-209. Herr, J. W., Chen, C. W., Goldstein, R. A., and Brogdon, J. N. "A Tool for Sediment TMDL Development on Oostanaula Creek." Proc., Total Maximum Daily Load (TMDL) Environmental Regulations: Proceedings of the March 11-13, 2002 Conference. Herr, J. W., Chen, C. W., Goldstein, R. A., Herd, R., and Brown, J. M. (2003). "Modeling acid mine drainage on a watershed scale for TMDL calculations." JAWRA Journal of the American Water Resources Association, 39(2), 289-300. Herr, J. W., Vijayaraghavan, K., and Knipping, E. (2010). "Comparison of Measured and MM5 Modeled Meteorology Data for Simulating Flow in a Mountain Watershed." Journal of the American Water Resources Association, 46(6), 1255-1263. Huber, W. C., Dickinson, R. E., and Barnwell Jr, T. O. (1988). "Storm water management model; version 4." Environmental Protection Agency, United States. Ji, Z. G., Morton, M. R., and Hamrick, J. H. "Modeling hydrodynamic and water quality processes in a reservoir." ASCE, 38. Ji, Z. G., Morton, M. R., and Hamrick, J. M. (2001). "Wetting and drying simulation of estuarine processes." Estuar Coast Shelf S, 53(5), 683-700. Jin, K. R., Ji, Z. G., and Hamrick, J. H. (2002). "Modeling winter circulation in Lake Okeechobee, Florida." Journal of Waterway, Port, Coastal, and Ocean Engineering, 128(3), 114-125. 124 Keller, A. (2000). "Peer Review of the Watershed Analysis Risk Management Framework (WARMF)?An evaluation of WARMF for TMDL applications by independent experts using USEPA guidelines." EPRI, Palo Alto, CA., Report, 1000252. Keller, A. (2001). "Peer Review of the Acid Mine Drainage Module of the Watershed Analysis Risk Management Framework (WARMF)?An evaluation of WARMF/AMD using USEPA guidelines." Electric Power Research Institute, Palo Alto, CA, Technical Report. Krause, P., Boyle, D. P., and Base, F. (2005). "Comparison of different efficiency criteria for hydrological model assessment." Advances in Geosciences, 5(89-97). Legates, D. R., and McCabe, G. J. (1999). "Evaluating the use of ?goodness-of-fit? measures in hydrologic and hydroclimatic model validation." Water Resources Research, 35(1), 233-241. Luo, D., Zhang, H., Zhang, F., Li, D., and Wang, X. "Applying WARMF model for SS spatial distribution in Chao Lake Basin, China: For example Hangbu-Fengle River watershed." Proc., 5th International Conference on Bioinformatics and Biomedical Engineering, iCBBE 2011, May 10, 2011 - May 12, 2011, IEEE Computer Society, IEEE Engineering in Medicine and Biology Society; Wuhan University; Fuzhou University; Nankai University; Overs. Chin. Sch. Environ. Prot. Assoc. (OCSEPA). Moriasi, D. N., Arnold, J. G., Van Liew, M. W., Bingner, R. L., Harmel, R. D., and Veith, T. L. (2007). "Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations." Transactions of the ASABE, 50(3), 885-900. Nash, J. E., and Sutcliffe, J. V. (1970). "River flow forecasting through conceptual models part I - a discussion of principles." Journal of Hydrology, 10(3), 282-290. 125 NC DENR (2009). "Falls Lake Watershed Analysis Risk Management Framework (WARMF) Development: Final Report." N. C. Department of Environment and Natural Resources, Raleigh, NC. Neilson, B., Horsburgh, J., Stevens, D., Matassa, M., Brogdon, J., and Spackman, A. (2003). "Comparison of Complex Watershed Models ?Predictive Capabilities: EPRI?s Watershed Analysis Risk Management Framework WARMF vs. USEPA?s Better Assessment Science Integrating Point and Nonpoint Sources BASINS/WinHSPF." Final Rep. Prepared for Utah Water Research Laboratory, Utah State Univ., Logan, Utah. Neitsch, S., Arnold, J., Kiniry, J., and Williams, J. (2005). "Soil and water assessment tool theoretical documentation, version 2005." Temple, TX. USDA Agricultural Research Service and Texas A&M Blackland Research Center. Parsons, J. E., Thomas, D. L., and Huffman, R. (2004). "Agricultural Non-Point Source Water Quality Models, their use and application." . (January 9, 2013). Rich, P. M., Weintraub, L. H., Chen, L., and Herr, J. "Climate Change Impacts on Hydrology and Water Management of the San Juan Basin." Proc., American Geophysical Union, Fall Meeting 2005, 1054. Riggs, T. L., Wilson, C. J., Ewers, M. E., Rich, P. M., and Weintraub, L. H. Z. (2005). "Decision Support for Water Planning: The ZeroNet Water-Energy Initiative." Impacts of Global Climate Change, 1-11. 126 Shrestha, S. (2010). "Simulations of Watershed Response to Land Use and Climate Change in the Saugahatchee Creek Watershed using the WARMF Model." M.S.C.E., Auburn University. Smakhtin, V. U. (2001). "Low flow hydrology: a review." Journal of Hydrology, 240(3-4), 147-186. Stringfellow, W., Herr, J., Litton, G., Brunell, M., Borglin, S., Hanlon, J., Chen, C., Graham, J., Burks, R., and Dahlgren, R. (2009). "Investigation of river eutrophication as part of a low dissolved oxygen total maximum daily load implementation." Water Science & Technology, 59(1), 9-14. Systech Engineering Inc (2005). "Creating a WARMF 6.1 Application Using a BASINS 3.1 Delineation: A User's Guide." USEPA (2004). "Better Assessment Science Integrating point and Nonpoint Sources (BASINS) Version 3.1: User's Manual." United States Environmental Protection Agency. USEPA (2007). "Model Report for Hangman (Latah) Creek TMDL Model Project." T. C. Group, ed., U.S. EPA. Vogel, R. M., and Fennessey, N. M. (1995). "Flow Duration Curves .2. A Review of Applications in Water-Resources Planning." Water Resour. Bull., 31(6), 1029-1039. Wang, P., Shariq, L., Montague, L., Kwaan, R., and Kella, V. (2004). "Developing a Nutrient Management Plan for the Napa River Watershed." Donald Bren School of Environmental Science and Management, University of California, Santa Barbara, CA. 127 Weintraub, L. H. Z., Chen, C. W., Goldstein, R. A., and Siegrist, R. L. "WARMF: A watershed modeling tool for onsite wastewater systems." Proc., 10th National Symposium on Individual and Small Community Sewage Systems, March 21, 2005 - March 24, 2005, American Society of Agricultural and Biological Engineers, 636- 646. 128