Copyright by Shiva Niazi 2000 Sensitivity Analysis of Time-Step in Modeling River and Aquifer Interaction by Shiva Niazi, B.S. Thesis Presented to the Faculty of the Graduate School of The University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of Master of Science in Engineering The University of Texas at Austin August 2000 Sensitivity Analysis of Time-Step in Modeling River and Aquifer Interaction Approved by Supervising Committee: Randall J. Charbeneau David R. Maidment to Mom and Dad v Acknowledgements I would like to express my gratitude to the many people who helped me during my time at the University of Texas at Austin. First, I would like to thank my advisors, Dr. Maidment and Dr. Charbeneau, for their guidance and encouragement during these past two years. I would also like to thank Dr. Maidment's research group, as well as Dan Opdyke, for patiently answering my technical questions and sharing their time with me. To the EWRE gang (Julie, Dave, Misti, Steph, Craig, Ann and Christine), you will be my favorite memory of my time here - thank you for your generosity, spirit and sense of humor. A heartfelt thanks goes to John, Shilpa and Ben for providing me a family away from home. I also would like to show my appreciation for Dr. Lawler and those in the Water Treatment in Developing Countries class (Katherine, Pablo, Caroline, Todd, Terri, Mary, Dave F., Dave A. and Andrea) for reminding me of why I wanted to become an engineer. Most importantly, I would like to thank my family for their support and for inspiring me to pursue higher education. August 10, 2000 vi Sensitivity Analysis of Time-Step in Modeling Aquifer and River Interaction Shiva Niazi, M.S.E. The University of Texas at Austin, 2000 Co-Supervisors: Randall J. Charbeneau and David R. Maidment ABSTRACT In modeling groundwater and surface water systems simultaneously, the issue of time-step length is important because of the difference in residence times of water in rivers and aquifers. To determine the effect of time-step length in modeling river and aquifer systems, a MODFLOW groundwater model of Milam, Lee and Bastrop counties Texas was dynamically linked to a surface flow routing model of the Colorado River. In the dynamic link between separate surface water and groundwater models, the output of one model is used to update the input of the other model in a cyclic fashion. Time-step length is defined as the length of time each model is allowed to run before updating the other model. A series of 32?day flood wave simulations was performed to determine the effect of averaging a highly fluctuating river discharge over 2, 4, 8 and 16-day time-steps. The results of this study suggest that time-step affects the quantity of water that the model predicts is exchanged between the river and aquifer. vii Table of Contents List of Tables........................................................................................................... x List of Figures..........................................................................................................xi Chapter 1: Introduction............................................................................................1 1.1 Background............................................................................................1 1.1.1 Observational Data .......................................................................2 1.1.2 Dynamically-Linking Models ......................................................2 1.1.3 Time Terminology........................................................................4 1.2 Research Objectives ..............................................................................6 1.3 Study Area .............................................................................................6 1.4 Summary................................................................................................7 Chapter 2: Literature Review ..................................................................................9 2.1 Introduction...........................................................................................9 2.2 Water Regulation...................................................................................9 2.2.1 Texas Water Law........................................................................10 2.3 Existing Models ...................................................................................12 2.3.1 MIKE-SHE .................................................................................13 2.3.2 MODFLOW-SURFACT............................................................15 2.3.3 Approximations of the Saint Venant's Equations .......................17 2.4 Time-step Analysis..............................................................................19 Chapter 3: MODFLOW Groundwater Model.......................................................21 3.1 Introduction.........................................................................................21 3.2 MODFLOW ........................................................................................22 3.2.1 Organization of MODFLOW Program ......................................23 3.3 River Package ......................................................................................27 3.4 Bureau of Economic Geology Model..................................................30 viii 3.4.1 Changes to the River Package Cells ...........................................33 3.4.2 Software Incompatibility............................................................35 Chapter 4: Surface Water Model...........................................................................37 4.1 Introduction.........................................................................................37 4.2 Observational Data ..............................................................................37 4.3 Manning's Equation.............................................................................44 4.4 Kinematic Wave Theory......................................................................46 4.5 Explicit Solution to the Kinematic Wave Equation............................50 4.5.1 Courant Condition......................................................................52 4.6 Determining Initial Flowrate in the Colorado River ...........................54 4.6.1 Data.............................................................................................54 4.6.2 GIS Processing of Data ..............................................................55 4.6.3 Upstream Flowrate .....................................................................60 4.6.4 Runoff and River Reach Lengths ...............................................62 4.6.5 Surface Water Model Calibration...............................................63 4.7 Excel-Based Surface Water Model .....................................................65 4.7.1 Surface Water Simulation...........................................................65 4.7.2 Surface Water Model..................................................................69 Chapter 5: Visual Basic Interface ..........................................................................76 5.1 Introduction.........................................................................................76 5.2 Overview .............................................................................................77 5.2.1 Input............................................................................................77 5.2.2 Dynamic Linkage .......................................................................79 5.2.3 Output .........................................................................................80 5.3 Structure ..............................................................................................81 5.4 Formatting of MODFLOW Files.........................................................84 5.5 Interface Enhancement ........................................................................84 5.5.1 Automation of MODFLOW Run ...............................................84 ix 5.5.2 Efficiency...................................................................................85 Chapter 6: Results..................................................................................................87 6.1 Introductio n .........................................................................................87 6.2 Spatial Variation in River Reaches......................................................88 6.3 Aquifer Head in Layer 1 ......................................................................91 6.3.1 Effect of Time-Step ....................................................................93 6.3.2 Effect of Distance from River ....................................................94 6.4 Aquifer Head in Layer 2 ......................................................................96 6.5 River Leakage Rates ............................................................................98 6.6 Approximation of Saint-Venant Equations .......................................105 Chapter 7: Conclusions ........................................................................................108 7.1 Introduction.......................................................................................108 7.2 Objectives ..........................................................................................108 7.3 Software.............................................................................................110 7.3.1 The Role of GIS........................................................................110 7.3.2 FORTRAN Compiler ...............................................................111 7.4 Preliminary Results ...........................................................................111 7.5 Recommendations .............................................................................112 Appendix A: Surface Water Model Data ...........................................................116 Appendix B: Interface Visual Basic Code..........................................................122 Appendix C: Simulation Results........................................................................136 References ..........................................................................................................177 Vita ....................................................................................................................178 x List of Tables Table 2.1 A Description of the MIKE-SHE Modules Excluding the Overland and River Module .........................................................................................14 Table 3.1 MODFLOW Packages and Their Purpose ............................................24 Table 4.1 GIS Data and Their Sources..................................................................55 Table 4.2 Data from USGS Gages Located Upstream of the Subject Area..........61 Table 5.1 The Nine Sub-Procedures of the "cmdOK_click" Procedure ...............82 Table 6.1 Volume of Groundwater Flow Entering in the Colorado River of the 32-Day Total Simulation Time ..................................................................103 Table 6.2 Volume of Water Leaked from River During the First 8, 16, 24 and 32 Days of Simulation ...............................................................................104 xi List of Figures Figure 1.1 Schematic of the Dynamically Linked Groundwater and Surface Water Models .................................................................................................3 Figure 1.2 Schematic of Model Timelines ..............................................................5 Figure 1.3 Location of Study Area..........................................................................7 Figure 3.1 Schematic of MODFLOW River Package ...........................................30 Figure 3.2 Layer 1 of the BEG Model and Revised BEG Model..........................32 Figure 3.3 Layer 1 of BEG Model in GMS Software ...........................................32 Figure 3.4 River Package Cells Used in the BEG Study.......................................34 Figure 3.5 River Package Cell Used in the Revised BEG Model.........................34 Figure 4.1 Location of the USGS Stream Gage and the Three Groundwater Wells .......................................................................................................................38 Figure 4.2 Location of USGS Gage 8159200 with Respect to Layer 1 of Revised BEG Model....................................................................................................39 Figure 4.3 Comparison of the Aquifer Head Deviation from the Average Head and the River Discharge at USGS Gage 8159200 .........................................40 Figure 4.4 Head Deviation in Well 584706 as a Function of River Discharge Deviation.......................................................................................................41 Figure 4.5 Comparison of Daily River Discharge at USGS Gage 8159200 and Aquifer Heads at Three Nearby Wells ..........................................................43 Figure 4.6 Schematic of the 8-Direction Pour Point Model..................................58 Figure 4.7 The Flow Direction Grid for the Subwatersheds Underlying Milam, Lee and Bastrop Counties..............................................................................59 Figure 4.8 Location of the Three USGS Gages Used in Determining the Upstream Flowrate in the Colorado River .....................................................................61 Figure 4.9 Finite Difference Nodes Along the Colorado River ............................63 Figure 4.10 Location of USGS Gage 8160400 with Respect to the Finite Difference Nodes...........................................................................................64 Figure 4.11 Time-Step Averaged River Discharge at the Inlet.............................67 Figure 4.12 Aquifer Cross Sections 1 and 2 in Layer 1 of the Revised BEG Model .......................................................................................................................69 Figure 4.13 A Section of the inlet2 Worksheet .....................................................72 Figure 4.14 A Section of the Par Worksheet........................................................73 Figure 4.15 A Section of the Q Worksheet ...........................................................75 Figure 5.1 Schematic of the Three Main Functions of the Interface.....................77 Figure 5.2 "frmInput" Form...................................................................................79 Figure 5.3 Flowchart of the "cmdOK_click" Procedure?????????... 83 Figure 6.1 River Discharge Along the Colorado River During Three 2-Day Time- Steps ..............................................................................................................89 xii Figure 6.2 River Discharge as Function of Simulation Time at Three Locations Along the River. ............................................................................................91 Figure 6.3 Location of Selected Cells in Cross Section 1 .....................................92 Figure 6.4 3-Dimensional View of the Selected Cells in Each Layer of Cross Section 1 ........................................................................................................92 Figure 6.5 Changes in Aquifer Head in the Cell Containing Reach 14 ................94 Figure 6.6 Aquifer Heads in Layer 1 of Cross Section 1 Using a 2-Day Time-Step .......................................................................................................................95 Figure 6.7 Aquifer Heads in Layer 2 of Cross Section 1 Using a 2-Day Time-Step .......................................................................................................................97 Figure 6.8 River Leakage Rate in Reach 14 Using Three Different Time-Steps..99 Figure 6.9 Sum of River Leakage Rates for All 35 Reaches Using Three Different Time-steps .....................................................................................................99 Figure 6.10 Aquifer Heads in Layer 1 of Cross Section 1 ..................................101 1 CHAPTER 1: INTRODUCTION 1.1 Background The successful management of water resources involves managing the two main components of a region?s water resources, namely groundwater and surface water. Surface water and groundwater are often managed separately (Lusk 1998), the fact that they are known to exchange water creates a strong incentive for the conjunctive management of these two resources. However, before the conjunctive management of surface water and groundwater can occur, it is imperative to determine how these two systems interact. One of the challenges in understanding the interaction of surface and ground water systems lies in their different time scales. Rivers, as a subset of surface water systems, have a much shorter residence time for water than do aquifers. Aquifers have much slower flow velocities and consequently may show slower changes in hydraulic head over time. Therefore, questions arise regarding the time-scale with which river and aquifers interact. Do the relatively fleeting river discharge fluctuations make an impact on the aquifer heads? If so, how are these changes distributed in time and space within the aquifer? 1.1.1 OBSERVATIONAL DATA To determine whether observational data could provide insight into comprehending a river and aquifer's effect on one another, a study was performed using data provided by a USGS stream gage and three nearby groundwater wells. 2 The wells and stream gage are located in Bastrop County, inside the study area shown in Figure 1.3. The study demonstrated the difficulty in using sparsely gathered observational data to determine the dynamic nature of aquifer and river interaction. Only a minimal level of correlation was observed between the aquifer and river hydraulic heads as a function of time. Furthermore, the spatial distribution of the river's impact on the aquifer was difficult to examine due to a lack of groundwater data. The details of this study are presented in Section 4.2, however from this analysis it was learned that a physically distributed model may provide a better understanding of an aquifer and river's interaction, as compared to relying on observational data alone. 1.1.2 DYNAMICALLY-LINKING MODELS In modeling the river and aquifer flow systems, the inherent difference in their time scale makes choosing an appropriate time-step difficult. A long time- step, which would be appropriate for modeling a groundwater system alone, might cause a loss of accuracy by over-averaging the river stage values. A short time-step, while good for modeling river systems, would substantially increase the computation time, and render the process inefficient for projecting water availability in the distant future. Therefore, it is hypothesized that an optimum time-step exists which would balance the need for accuracy and an efficient modeling system. In this project, a dynamic link was created between a surface water and groundwater model to help assess the role of time-step in the modeling of the two systems. In the dynamic link between separate surface water and groundwater 3 models, the output of one model is used to update the input of the other model in a cyclic fashion, as shown in Figure 1.1. A physically distributed model of a river was created specifically for this study in Microsoft Excel. A calibrated MODFLOW groundwater model that was developed by the University of Texas Bureau of Economic Geology was used in this research. As seen in Figure 1.1, the parameters that are exchanged between the two models are the river hydraulic head, HR and the lateral aquifer inflow into the river, qa. Figure 1.1 Schematic of the Dynamically Linked Groundwater and Surface Water Models 1.1.3 TIME TERMINOLOGY Figure 1.2 shows a schematic of the time terminology that is used in this thesis. For the purpose of this study, the time-step of interest is the Visual Basic Interface time-step shown in Figure 1.2. In a dynamic linkage of two models, the interface time-step length corresponds to the length of time each model is allowed 4 to run before updating the other model. Therefore, the run-time of a model (i.e. the length of time each model is run) is equivalent to the interface time-step. However, there are also separate time-steps intrinsic to the surface water and groundwater models. When discussing these time-steps the words groundwater or surface water will precede time-step to indicate its limited applicability to the scope of that one model. As shown in Figure 1.2, the surface water model consists of multiple surface water time-steps. However, due to the implicit finite difference scheme used by MODFLOW, the relatively short run-time of the MODFLOW model consists of only one groundwater time-step. Typical MODFLOW models that are used to model groundwater movement alone, however, often require multiple time-steps. The total simulation time is defined as length of time for which the dynamically linked models have been run. Therefore, the total simulation time is the summation of all of the interface time- steps. 5 Figure 1.2 Schematic of Model Timelines 1.2 Research Objectives There are three primary objectives for this research. The first is to build a simple but appropriate surface water model that may be used to simulate flow in real rivers, such as Colorado River. This model needs to be easily manipulated by the interface that is facilitating the dynamic linkage. The second objective is to design the code for the interface that dynamically links the models. In this project, the interface was written in the Visual Basic Language. The interface needs to be able to run each model in an alternating fashion and transfer the 6 necessary information between the models. Lastly, the effect of time-step in modeling the interaction between the groundwater and surface water systems is analyzed. To do so, a highly time-varying river discharge caused by a flood wave is routed through the domain using the kinematic wave equation. The river discharge is averaged over varying time-steps to determine if a relationship exists between the interface time-step and the river and aquifer interaction. 1.3 Study Area The study area chosen is located in Milam, Lee and Bastrop counties in Texas and is shown in Figure 1.3. The two main rivers in this area are the Colorado River in the west and the Brazos River in the east. The main aquifer analyzed in this region is the Carrizo-Wilcox aquifer. This area has recently undergone a groundwater availability investigation by the Bureau of Economic Geology (Dutton 1999). The groundwater model was therefore available from the Bureau of Economic Geology for use in this study. 7 Figure 1.3 Location of Study Area 1.4 Summary This research attempts to determine the interface time-step appropriate for modeling groundwater and surface water interaction. The case study considers a specific site in Texas. In this research a groundwater and surface water model are dynamically linked. In Chapter 2, the literature review for this research is presented. The groundwater and surface water models are described in Chapters 3 and 4, respectively. The description of the Visual Basic interface is left to Chapter 5. Lastly, in Chapter 6, the results of the different simulations with 8 varying time-step are reported and conclusions and recommendations are addressed in Chapter 7. 9 CHAPTER 2: LITERATURE REVIEW 2.1 Introduction This chapter discusses the available literature in the field of surface and subsurface water modeling. In Section 2.2 the current state of water regulations, as it pertains to both surface and subsurface systems in Texas, is examined. Next, Section 2.3 explains the different types of models used for water management purposes, with a detailed look at the MIKE-SHE and MODFLOW-SURFACT models. Lastly, previous parameter sensitivity analyses are examined and the possible contributions of this research are considered. 2.2 Water Regulation As water demand continues to increase with population growth, water management practices that promote water conservation become increasingly necessary. Unlike the Clean Water Act that sets federal limits on the contamination of all ?navigable? waters in the U.S., no such federal law exists for the withdrawal of water (Vance 1996). Instead, individual States are given the right to protect their own water resources as they see fit. This lack of federal regulations creates a varied approach to State resource management practices, with the result that some States are better at passing laws to conserve their resources than others (Vance 1996). 10 2.2.1 TEXAS WATER LAW As compared with other state codes, Texas Water Law has been slow to keep up with more advanced water management practices. Before the passing of Senate Bill 1, the extensive water management bill, in 1997, Texas was one of three western states without a state drought plan (Hubert 1999). Furthermore, Texas is the only western state that still abides by the rule of capture with respect to groundwater pumping (Vance 1996). The rule of capture states that landowners have property rights to the water below their land, and therefore can withdraw groundwater at their discretion1 (Lusk 1998). As a result, there is no motivation for landowners to conserve groundwater under this rule, leading to the overuse of groundwater resources (Lusk 1998). Unlike groundwater, surface water withdrawals are governed by the State of Texas under the system of prior appropriation. Under prior appropriation, priority is given to surface water rights based on the dates the permits were issued under the doctrine of ?first in time...first in right.? (Lusk 1998) In a conflict between two water rights during a water shortage, the senior water right, namely the one that was issued first, can exercise its full water rights before the junior water right can use the water. Under this system, the power to decrease water usage during droughts is left within the hands of the State (Lusk 1998). Paradoxically, Texas does not regulate the pumping of groundwater wells that could be robbing baseflow to a nearby river. Pumping of such groundwater would lead to a decrease in the river?s flowrate, as would a surface water 1 In Texas, the rule of capture has been modified; groundwater pumping can be curtailed if it is proven 1.) to be malicious or 2.) to constitute willful waste. (Vance 1996) 11 withdrawal. However, under the rule of capture, this pumping could not be regulated by the State (Lusk, 1998). Therefore, this system of regulation can create inequities between groundwater and surface water users. With the passing of Senate Bill 1, Texas Water Law has taken large strides towards managing its scarce water resources (Hubert 1999). Although there is still no state regulation on groundwater pumping, Senate Bill 1 did give more power to local communities to alter the rule of capture, through revising its legislation toward groundwater conservation districts (Hubert 1999). Should they choose, groundwater conservation districts have the power to deny groundwater well permits based on numerous criteria, including if the ?proposed use of water unreasonably affects existing ?surface water resources? (Texas Water Code 36.113). Also due to Senate Bill 1 (Texas Water Code 16.012), an effort has begun to assess the extent of the state?s surface water and groundwater resources through the Water Availability Modeling (WAM) and the proposed MODFLOW-based Groundwater Availability Modeling (GAM) projects (Mace and Mullican 2000). A complete water resources management plan requires an understanding of how the groundwater and surface water systems affect one another. Texas water resource planning and management stands, therefore, to gain from models that can predict the dynamics of groundwater and surface water interaction. 12 2.3 Existing Models Since the 1980?s, various states have been making use of joint surface/subsurface modeling systems to help create legislation that would equitably manage their water resources (Sophocleous 1995; Mueller 1993). Models can provide information where data is not available, for example to project water availability in the future. Models typically fall under two categories: lumped conceptual and physically based distributed (Refsgaard 1997). Physically based distributed models can represent spatially varying parameters that are based on physical characteristics of the system. Lumped conceptual models treat complex physical systems, such as a watershed, as an integrated unit (Refsgaard and Knudsen 1996). The advantage of using a physically based distributed model is that localized changes, such as changes in land-use, can be modeled more effectively. Furthermore, in theory, physically based distributed models require less time-series data for calibration (Abbott 1986). In recent years, there has been an emphasis on linking physically distributed models to a Geographic Information Systems (GIS). GIS, as a powerful graphical tool, has helped in the visualization of model output (Orzol 1993). Furthermore, models that were previously created in arbitrary coordinates can now, with a GIS, be geographically referenced. This ability to view and manipulate disparate data sets within one integrated system increases the ease and efficiency with which physically based data can be gathered and formatted for input into the model (Hinaman 1993; Richards 1993). 13 Computer advances have also created increasingly sophisticated models that are able to integrate complicated physical processes such as runoff, river flow, subsurface flow, evapotranspiration and contaminant transport. Two such models will be discussed below: MIKE-SHE and MODFLOW-SURFACT. 2.3.1 MIKE-SHE MIKE-SHE is the most recent development of the Systeme Hydrologique Europeen hydrological model created in a joint project by the Danish Hydrological Institute, the British Institute of Hydrology and the French consulting company SORGEAH (Abbott 1986). Like MODFLOW, MIKE-SHE consists of separate modules that model different aspects of the hydrological cycle (Danish Hydraulic Institute 2000b). The basic modules include the following: Overland and River Module, Evapotranspiration Module, Unsaturated Zone Module, Saturated Zone Module and the Irrigation Module. Note that there are additional modules that can be purchased separately. The basic modules, excluding the Overland and River Module, which will be discussed in more detail in the following paragraph, are described in Table 2.1. 14 Table 2.1 A Description of the MIKE-SHE Modules Excluding the Overland and River Module MODULE DESCRIPTION Evapotranspiration (ET) Models rain interception and evapotranspiration by either: i. Rutter model (interception) and Penman- Monteith Equation (evapotranspiration) ii. Kristensen-Jensen model Unsaturated Zone Module (UZ) Models vertical flow in unsaturated zone of surface by either: i. Richard?s Equation ii. Gravity flow Saturate Zone (SZ)2 Models saturated flow in subsurface. Solver methods: i. Pre-conditioned Conjugate Gradients (PCG) ii. Modified Gauss-Seidel Irrigation (IR) Models irrigation management: Highly flexible, can place priority on water source The Overland and River Module models the overland runoff and river flows in tandem. The two-dimensional diffusive wave approximation of the Saint-Venant equations and Manning equation are used for modeling the overland flow. The river can be modeled in two levels of complexity, using 1) a one- 2 It should be noted that the river/aquifer interaction is specified in the SZ Module. The river/aquifer interaction can be modeled assuming: 1.) the river is in full contact with the aquifer or 2.) a low permeability layer separates the river and aquifer. 15 dimensional diffusive wave approximation for the Saint-Venant equations or 2) MIKE-11, a one-dimensional river hydraulic model distributed by the Danish Hydraulic Institute. MIKE-11 is also structured in a modular fashion, with modules that describe river flow (Hydrodynamic Module), contaminant transport (Advection-Dispersion Module) and biological processes (Water Quality Module). The Hydrodynamic Module solves the Saint-Venant equations for open channel flow for fully dynamic, diffusive, kinematic and quasi-steady state waves. Furthermore, this module can model flows around a variety of structures including broadcrested weirs and culverts (Danish Hydraulic Institute 2000a). 2.3.2 MODFLOW-SURFACT The MODFLOW-SURFACT modeling system is a product of HydroGeologic, Inc. MODFLOW is a modular 3-Dimensional finite difference modeling system of saturated subsurface flow that was created by the United States Geologic Survey (USGS). MODFLOW-SURFACT integrates an enhanced version of MODFLOW-96 with packages that model overland flow, channel flow and contaminant transport. Open?channel flow and overland flow are modeled using the 1-D Channel Flow (CHF1) and 2-D Areal Overland Flow (OLF1) Packages, respectively. The OLF1 Package uses the two-dimensional diffusive wave approximation to model overland flow. The package provides an extra layer of nodes that are located above the aquifer layers that are modeled by MODFLOW. These overland flow nodes are able to exchange water with the first active layer of groundwater nodes directly below them via leakage through the soil surface. The CHF1 Package uses the one-dimensional diffusive wave 16 approximation to model flow in an open channel. Two types of channel cross sections are supported in the CHF1 Package: a wide rectangular channel and a trapezoidal channel. The aquifer and river interaction is modeled using Darcy?s Law, which is described in more detail in Chapter 3 of this thesis (HydroGeologic 1999). Both packages provide the following five types of boundary conditions: ? First type ? Areal recharge ? Sources and Sinks ? Evaporation ? Zero-depth gradient and critical depth gradient The following is a brief explanation of the five boundary conditions listed above. The First Type boundary condition is equivalent to the constant head boundary option that is offered by MODFLOW, in which the heads at the boundary are kept at the initial head value throughout the simulation. The second boundary condition (Areal Recharge) applies a recharge rate to a horizontal area in the OLF1 Package and the channel surface area in the CHF1 Package. This boundary condition contains a maximum depth constraint that can limit the recharge rate. The Sources and Sinks boundary condition and the Evaporation boundary condition provide net fluxes and an areal sink to the boundary node, respectively. The Sink and Evaporation conditions are subject to a non-negative depth constraint and, therefore, cannot cause the hydraulic head at the node to 17 drop below the bed elevation. The last boundary condition, the Zero-Depth Gradient and Critical Depth Gradient, simulates the condition at the bottom of a hill or at the downstream end of a river reach. The Zero-Depth Gradient condition causes the slope of the water surface to equal the riverbed slope, while the Critical Depth Gradient condition forces the water depth at the boundary to be equivalent to the critical depth. Both the overland and channel flow algorithms are fully integrated into an implicit system of matrix equations and are solved during each time-step. The overland and open channel nodes modeled separately by the two packages, exchange water via equations for a free-flowing weir and submerged weir. The free-flowing weir equation is used under the condition that the hydraulic head within the channel is lower than the elevation of the channel bank. Alternatively, a submerged weir condition occurs when the hydraulic head within the channel is higher than the channel bank elevation (HydroGeologic 1999). Other features of the enhanced MODFLOW modeling system include a 3- dimensional vadose zone transport addition to the Block-Centered-Flow Package and an advanced time-stepping mechanism in the ATO4 package. The advanced time-stepping package can modify the groundwater model time-step depending on the computer?s available memory so as to have the solver algorithm more easily converge to a solution (HydroGeologic 1999). 2.3.3 APPROXIMATIONS OF SAINT VENANT'S EQUATIONS Although both the MIKE-SHE and MODFLOW-SURFACT modeling systems use the diffusive wave approximation of the Saint-Venant equations, 18 Lighthill and Whitham (1955) showed that the main part of a natural flood wave approximates a kinematic wave (a simplified version of the diffusive wave). Because of the flood wave simulations investigated in this study, the kinematic wave approximation is believed to be suitable. However, in instances where changes in river discharge are based solely on lateral inflow and not influenced by a flood wave, Vieira (1983) investigated the accuracy of different approximations of the Saint-Venant equations based on two parameters: the kinematic wave number, k, and the Froude number, Fo. Equations 2.3.1 and 2.3.2 define the two parameters. The approximations were solved under two lower boundary conditions: the zero-depth-gradient and critical-flow. These two conditions are identical to the zero-depth gradient and critical-depth gradient boundary conditions used in the MODFLOW-SURFACT modeling system. In the study, Vieira compared the implicit finite difference solutions to the kinematic, diffusive and gravity wave equations to those of the full Saint-Venant equations. The results of the study showed that for values of k much greater than 50, the kinematic wave approximation can be used with either boundary condition. However, for k values between 5 and 20, the kinematic and diffusive wave approximations are applicable depending on the Froude number. For example, a flow regime with the parameters k = 20 and Fo < 0.5 would require the diffusive wave approximation as opposed to the more simplistic kinematic wave approximation. The boundary conditions become significant only when k < 5. 19 3/1 24 3 sin ??? ? ??? ?= qC Lgk q (2.3.1) 2/1tan ??? ? ??? ?= gCFo q (2.3.2) where, k= kinematic wave number Fo = Froude number g = gravitational acceleration [L2/T] L = length of river reach [L] ? = constant angle of the slope R = hydraulic radius [L] q = lateral runoff [L/T] 2.4 Time-Step Analysis The hydrologic modeling systems discussed above are powerful in their ability to model complex physical processes, such as river and aquifer interaction, simultaneously. However, to use these modeling systems effectively it is important to choose spatial and time variables that model the system of interest to the desired degree of accuracy. There have been studies done on the changes that various physical and structural parameters of the coupled surface water and groundwater models can make on the river and aquifer hydraulic heads 20 (Refsgaard 1996; Bathurst 1986; Perkins 1999). However, attempts to look at the effect of time-step have been, at best, brief. This research, therefore, aims to continue the work of previous sensitivity analyses of surface and groundwater models with an emphasis on examining the effects of changing the time-step. Furthermore, the linking of the models is done in a GIS context. In this way, the benefits that GIS can provide in linking two separate physically distributed models are examined. 21 CHAPTER 3: MODFLOW GROUNDWATER MODEL 3.1 Introduction The groundwater model, MODFLOW, was chosen to model the groundwater system in this project for two reasons. Firstly, MODFLOW is well established and widely used in the United States in the fields of civil and geotechnical engineering, and in hydrogeology, making this research directly relevant to the work of many industries modeling subsurface flow. Secondly, the model contains a module known as the River Package that is able to model the interaction between the river and aquifer, albeit in a non-dynamic fashion. Therefore, MODFLOW has the means to dynamically link the surface and subsurface models embedded in its program. The following section of this chapter describes the general capabilities of MODFLOW and its organization. In Section 3.3, the River Package and the underlying theory of its code are described in more detail. Finally, the specific parameters that describe the University of Texas Bureau of Economic Geology (BEG) MODFLOW model of the area within the Milam, Lee and Bastrop counties are discussed in Section 3.4. 3.2 MODFLOW The Three-Dimensional Modular Ground-Water Flow Model (MODFLOW) was created by the United States Geological Survey (USGS) in 22 1983 and since then has been updated numerous times, the latest version being MODFLOW-96. Since its inception, MODFLOW has become one of the most frequently used groundwater models in both academia and industry. MODFLOW is written in FORTRAN and uses a block-centered finite difference technique to solve the mass conservation equation that describes subsurface flow (Equation 3.2.1). In many "real-life" systems, complexities such as irregular model geometry, heterogeneous parameters, complex boundary conditions or some combination thereof, often make analytical solutions impossible. In such cases, MODFLOW and other models can provide numerical solutions. MODFLOW solves Equation 3.2.1 by employing the finite-difference method in a time- iterative fashion. In the most general sense, MODFLOW determines the aquifer hydraulic heads as a function of time based on boundary and initial conditions and stresses on the aquifer being modeled. Aquifer stresses include well pumping, interaction with rivers, and area recharge as well as others. ')()()( WzhKzyhKyxhKxthS zzyyxxS +????+????+????=?? (3.2.1) where, Ss = specific storage [L-1] Kii = principal components of the hydraulic conductivity tensor [L/T] h = hydraulic aquifer head [L] W?= source strength [Volume/(VolumeT)] 23 3.2.1 ORGANIZATION OF MODLFOW PROGRAM The MODFLOW model is organized into a main program and several independent modules called packages. Some of the packages are described in Table 3.1. The modular organization of the model allows the user to choose the packages that are needed to describe the system being modeled. For instance, depending on whether wells are located in the domain, the Well Package can be turned on or off. Therefore, unnecessary packages are ignored and do not increase the run time of the model. 24 Table 3.1 MODFLOW Packages and Their Purpose (Charbeneau 2000) Package Purpose Basic Package Handles tasks that are required for each simulation. Including specification of boundaries, determination of time-step length, establishment of initial conditions, and printing of results. Block-Centered Flow Package Calculates hydraulic conductance and external source terms of finite-difference equations that represent flow from cell to cell and storage. River Package Stress package. Adds terms representing flow to rivers to the finite-difference equations. Recharge Package Stress package. Adds terms representing diffuse recharge to the finite-difference equations. Well Package Stress package. Adds terms representing flow to wells to the finite-difference equations. Drain Package Stress package. Adds terms representing flow to drains to the finite-difference equations. Evapotranspiration Package Stress package. Adds terms representing evapotranspiration to the finite difference equations. General-Head Boundary Package Stress package. Adds terms representing general- head boundaries to the finite-difference equations. Solution Procedure Package MODFLOW (1996) supports preconditioned- conjugate gradient, strongly-implicit, slice- successive over relaxation, and direct solver using diagonal ordering procedures. In MODFLOW-96, the main program receives the Name file or the file consisting of the names and unit numbers of the different packages being used for 25 the simulation. Table 3.1 describes the different contents of the packages in MODFLOW-96. The most basic model (i.e. with no aquifer stresses) needs a minimum of the Basic, Block-Centered-Flow (BCF), Solution Procedure, and Output Control Packages to be defined by the user. The following paragraphs provide a brief description of the different packages in MODFLOW required for a basic model. In the Basic Package, the number of rows, columns and layers are specified. The initial hydraulic heads as well as the groundwater stress period and time-step are also specified in the Basic Package. Stress periods coincide with periods where parameters specifying the aquifer stresses (such as pumping rate, river stage, etc.) and boundary conditions must be held constant. Stress periods are further broken down into groundwater time-steps. The finite difference equations are solved iteratively for each of these groundwater time-steps. Typical MODFLOW groundwater stress periods can range from months to years, depending on the objectives of the model. As discussed before, the duration of groundwater time-steps are very much dependent on the solver algorithm and are therefore difficult to generalize. The BCF Package contains parameters, such as the hydraulic conductivity, aquifer types and row and column spacing, which describe the cell-to-cell flow. The spacing does not need to be uniform across the grid, and often to decrease the computation time, a modeler will have large grid cell sizes near the boundaries of the study area where accuracy is less important. 26 Thirdly, the Solution Procedure Package contains the information on what kind of solution method will be applied to solve the finite difference equations. Examples of two different solution methods are the Strongly Implicit Procedure (SIP) and Slice-Successive Over Relaxation (SSOR) algorithms. Certain solution methods could fail to solve the finite difference equations because the model will be unable to converge to a solution. Convergence, however, is not only dependent on the solution method but also relies on other parameters such as MODFLOW's time-step length. Therefore, by having different solution methods, the modeler is given more freedom in choosing values for parameters, such as the groundwater time-step length and various iteration parameters, which also influence the possibility of a solution. Finally, the Output Control Package specifies the format and content of the output files. There are two formats for the output of a MODFLOW simulation: a text file and a binary file. The text file, known as the listing file, lists each computational process performed by MODFLOW as it occurs during the simulation. The listing file also may contain, among other things, the aquifer heads and drawdowns at any groundwater time-step specified by the user. Through the Output Control Package, the user can also save the river leakage rates for each River Package cell to the listing file. This last capability is important in facilitating the dynamic linkage of the MODFLOW and surface water models. Unlike text files, binary files can?t be opened in word processing software, and require a GUI to display their contents. Common binary file outputs include 27 the aquifer head file (*.hed), the drawdown file (*.drn) and the cell-to-cell flow file (*.ccf). The cell-to-cell flow contains a water budget for each grid cell by recording the amount of water that flows through a cell during a specified groundwater time-step. MODFLOW uses binary files as a means of saving output as well as a means of entering input. Therefore, the output of one MODFLOW simulation can be used as the input for another simulation via binary files. For this research, aquifer heads contained in the *.hed file resulting from a MODFLOW run in one interface time-step are used as the starting aquifer heads in the following interface time-step. In this way, the *.hed file creates continuity between the disparate MODFLOW runs in each interface time-step. To summarize, the modeler can designate whether a specific output, such as aquifer heads, is ?printed? to a text file and/or ?saved? to a binary file by changing the contents of the Output Control Package. 3.3 River Package MODFLOW consists of a River Package that models the water influx into or drainage out of the aquifer from overlying rivers. Figure 3.1 depicts how MODFLOW models the river-aquifer interaction. This interaction is based on Darcy's Law where the flowrate of water between the river and aquifer is directly proportional to the hydraulic head difference between the two. The exact form of Darcy's Law used by MODFLOW, which describes the river/aquifer interaction is: 28 Wi,j,k = Ci,j,k (HRIV - hi,j,k) for hi,j,k > HBOT (3.3.1) Wi,j,k = Ci,j,k (HRIV -HBOT) for hi,j,k < HBOT where, Wi,j,k = aquifer recharge rate [L3/T] Ci,j,k = riverbed conductance [L2/T] HRIV = hydraulic head in river [L] HBOT = elevation of the riverbed bottom [L] hi,j,k = hydraulic head in aquifer [L] i,j,k refers to the parameter value in row i, column j and layer k Note that the Darcy Law equation is defined with respect to two conditions. Under the first condition, the hydraulic head in the aquifer is greater than the bottom of the riverbed, causing the flowrate to be partially controlled by the head in the aquifer. However, once the hydraulic head in the aquifer falls below the riverbed elevation, the aquifer recharge rate is independent of the hydraulic head in the aquifer. Under this second condition, the recharge rate is constant and at its maximum value. In this research, the second condition was never met. The conductance term, Ci,j,k, is a function of the physical parameters of the river and is defined in Equation (3.3.2). In addition, Figure 3.1 shows a schematic describing the river/aquifer interaction as modeled by the River Package. 29 M wLKC kji =,, (3.3.2) where, w = width of river [L] L = length of river in cell i,j,k [L] M = riverbed thickness [L] K = hydraulic conductivity of the riverbed [L/T] In general, conductance values are very difficult to measure or assess with certainty. Although topographic maps or GIS coverages can lend some assistance in assigning the length and width of a river, the two riverbed parameters are difficult to ascertain. Often, to attain some estimate of the hydraulic conductivity, assumptions must be made about the texture of the riverbed material. Hydraulic conductivities of clays and silt soils are sometimes used for lack of better information. Overall, little riverbed data exist because soil samples are not typically taken from the riverbed for the purpose of determining the riverbed conductivity. 30 Figure 3.1 Schematic of MODFLOW River Package 3.4 Bureau of Economic Geology Model The MODFLOW groundwater model that was used in this project was developed by the University of Texas Bureau of Economic Geology (BEG). The objective of the BEG project was to determine the potential hydrologic impact of the purchase of groundwater from Milam, Lee and Bastrop counties by San Antonio on wells of private users in the vicinity. The groundwater model provided by the BEG is referred to as the "BEG model," while the altered version of the BEG model used in this study will be referred to as the ?revised BEG model.? This section will describe only the BEG model parameters that are relevant to the scope of this study. A more complete description of the BEG model is presented by Dutton (1999). 31 Like the BEG model, the revised BEG model consists of one confined aquifer layer (Layer 1) and four confined/unconfined aquifer layers. A confined/unconfined layer includes both the outcrop and downdip areas of an aquifer. Layer 1 of the groundwater model area is shown in Figure 3.2. The River Package module was used to describe the influence of the Colorado, Brazos and Yegua Creek Rivers on the Carrizo-Wilcox aquifer. Figure 3.2 delineates the grid cells that were designated as River Package cells in Layer 1 of the BEG model. As discussed in the following subsection, the revised BEG model altered the location of some of the Colorado River Package cells. The remaining four layers of the groundwater model also contain River Package cells. These cells occur where one of the three rivers cross the outcrop area of the lower layer. These River Package cells were neither altered nor included in the dynamic linkage. 32 Figure 3.2 Layer 1 of the BEG (and revised BEG) model. The turquoise and light blue areas signify the Carrizo-Wilcox aquifer's outcrop and downdip areas, respectively. The Colorado and Brazos Rivers are shown in dark blue and green. The bright green lines coincide with the boundaries of Milam, Lee and Bastrop counties. Figure 3.3 Layer 1, as used in the BEG study, seen in the GMS software. The blue crosses indicate the River Package cells. The constant head cells are displayed with orange dots. 33 3.4.1 CHANGES TO THE RIVER PACKAGE CELLS For the purpose of this study, only the grid cells underlying the Colorado River in Layer 1 were linked to the surface water model. Furthermore, these cells were rearranged from locations shown in Figure 3.4 to those in Figure 3.5. As can be seen from Figure 3.4, the Rf1 file of the Colorado River did not directly overlie many of the original River Package cells, making the river length in each cell impossible to assess. Even after this alteration, two cells in a winding section of the Colorado River posed hydrologic problems when modeled separately. Because each River Package cell in MODFLOW can have only one river head value, the two separate river reaches that cross grid cell 32 were forced to have identical head values. Therefore, for the sake of continuity, the river segment in cell 31, which joins the two segments in cell 32, also has the same river head. 34 Figure 3.4 (top) shows the River Package cells representing the Colorado River that were used in the BEG study. The cells were modified to those shown in Figure 3.5 (bottom) to better fit the Colorado River. The red lines indicate cells 31 and 32 that were modeled as one reach in the surface water model. To stay consistent with the BEG model, the BEG riverbed elevations were primarily adopted for this study. Cells that were designated as River Package cells in both the BEG and revised BEG models, kept their assigned BEG riverbed elevations. However, cells that were newly designated as River Package cells in 35 the revised BEG model were assigned elevations that were interpolated from the BEG riverbed elevations. River conductance terms (CR) for the Colorado River Package cells were designated as 48,000 ft2/d (0.052 m2/s) by the BEG. To confirm that these values were still reasonable after the alterations to the River Package cells, the hydraulic conductivities of the riverbed were back - calculated using Equation 3.3.2. The river reach lengths were determined using a Geographic Information System, while the river width was taken from the BEG study to be 76.2 meters. Lastly, the riverbed thickness was assumed to be 10 cm. The resulting riverbed conductivities were found to range between 1E-06 and 1E-07 cm/s. This range of hydraulic conductivities is characteristic of silt or loess (Charbeneau 2000), a soil type that is often deposited on riverbeds. Therefore, the conductance terms of the BEG model were determined to be suitable for this study. 3.4.2 SOFTWARE INCOMPATIBILITY It should be noted that some packages supported in MODFLOW-based software such as Visual MODFLOW and Groundwater Modeling System (GMS) are not supported by MODFLOW-96 as provided by the USGS. The BEG model was created using Visual MODFLOW. However, to simplify the interaction between the Visual Basic interface and MODFLOW, the original MODFLOW-96 executable as downloaded from the USGS website was used for this study. Therefore, three packages, the Horizontal Flow Barrier Package (HFB), Constant Head (CH) and the WHB Solver Package created for the BEG model could not be 36 used during the simulations presented in this study. However, this loss in modeling capability is not expected to affect the results of this study. 37 CHAPTER 4: SURFACE WATER MODEL 4.1 Introduction For this study, a surface water model was created with the use of a Geographic Information System (GIS) and Microsoft Excel (Excel). The model is based on Manning's equation and kinematic wave theory, and serves to quantify the change in river stage as a function of the groundwater recharge rate. Although surface water models abound, a simple model in Excel has the advantage of being easily manipulated by the Visual Basic interface. As mentioned in Chapter 1, the possible relationship between observational river and groundwater data was investigated in a brief study. The study's results and limitations are presented in Section 4.2. Next, the assumptions and theory of both Manning's equation and kinematic wave theory are presented in Section 4.3 and 4.4, respectively. Section 4.5 describes the numerical solution to the kinematic wave equation used in this study. The methodology employed in determining the initial flowrate in the Colorado River will be discussed in Section 4.6. Lastly, Section 4.7 discusses the structure of the Excel surface water model. 4.2 Observational Data To determine whether observational data could provide insight into comprehending a river and aquifer?s effect on one another, a USGS gage measuring streamflow in the Colorado River at Bastrop and three aquifer 38 observation wells near the river were chosen for analysis. Figure 4.1 shows the relative locations of the USGS stream gage 8159200 and of the three groundwater wells (5854801, 5854706, 5862603) used in this study. The USGS stream gage is located directly upstream of the BEG MODFLOW modeling domain, as shown in Figure 4.2. The USGS streamflow data was downloaded from the USGS web site (http://waterdata.usgs.gov/nwis-w/TX/). Similarly, the hydraulic head values for the aquifer were provided by the Texas Water Development Board web page (http://www.twdb.state.tx.us/data/groundwater/groundwater_toc.htm). Figure 4.1 Location of the USGS Stream Gage and the Three Groundwater Wells 39 Figure 4.2 Location of USGS Gage 8159200 with Respect to Layer 1 of Revised BEG Model Aquifer hydraulic head and river discharge data were analyzed between the months of April, 1981 and June, 1985. This period was chosen because it contained data that had been recorded consistently at least once every two or three months. Figure 4.3 shows a graph of the aquifer hydraulic head deviation in the three wells and the river discharge at the stream gage as a function of time. In this study, the head deviation (?h) corresponds to the difference between the instantaneous aquifer head and the time averaged aquifer head. In Figure 4.3, the aquifer head deviation is plotted on the primary axis on the left, while the river?s flowrate as measured by the USGS gage, is plotted on the secondary axis on the right. Based on visual analysis, there appears to be some correlation between the river discharge and aquifer head deviation. 40 Comparison of Groundwater Head Deviation from Average and Flowrate at USGS Gage 8159200 -4 -3 -2 -1 0 1 2 3 4 5 Apr-81 Sep-81 Jul-82 Nov-82 Apr-83 Jun-84 Jun-85 Date Groundwater Head Deviation from Average (ft) 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 River Discharge (cfs) 5854706 5854801 5862603 flowrate Figure 4.3 Comparison of the Aquifer Head Deviation from the Average Head and the River Discharge at USGS Gage 8159200 To quantify the level of correlation between the aquifer and river fluctuations, the head deviation (?h) from the average was plotted against the river discharge deviation from the average (?Q). Similarly, the river deviation from average (?Q) is defined as the difference between the instantaneous river discharge and the time averaged river discharge. Figure 4.4 presents these values for well 5854706. A linear regression of this data produced an R2 value of 0.4088. A similar analysis of wells 5854801 and 5862603 yielded R2 values of 0.1564 and 0.0085, respectively. Overall, the R2 values seemed to indicate that 41 there was moderate to no correlation between the changes in river flowrate and aquifer head near the river. Dh in Well 5854706 vs. DQ at USGS Gage 80159200 y = 0.0003x - 0.1244 R2 = 0.4088 -4.00 -3.00 -2.00 -1.00 0.00 1.00 2.00 3.00 4.00 5.00 6.00 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 DQ (cfs) DDh (ft) Figure 4.4 Head Deviation in Well 584706 as a Function of River Discharge Deviation From Figure 4.3, there appears to be a similar pattern in hydraulic heads in all three wells in response to a flood wave occurring between the months of May and June of the year 1981. Figure 4.5 presents daily river discharges and aquifer heads during the flood wave period. As displayed in the figure, the hydraulic heads in all three wells increased in accordance with the higher river discharge rates. Not surprisingly, well 5854706, which has the highest correlation of the three wells, also displays the highest change in head of nearly 3 ft. Similarly, the 42 smallest increase in hydraulic head (approximately 0.5 ft.) occurs in well 5862603, which also has the lowest R2 value. Therefore, in response to the flood wave, the head change observed in the wells ranged from 0.5 ft. to 3 ft. over a period of one month. However, due to the lack of data, the lateral distance from the river and the well screening intervals do not appear to provide insight to the spatial distribution of the river and aquifer interaction. Well 5854706 has the highest R2 value of the three wells and is the closest well to the river, however its screening interval is the lowest of the three wells at 400-440 ft. below the top of the casing. Well 5854801 shows a medium range of correlation, however it is the farthest from the river. The screening interval for this well is in the middle range at a depth of 300-360 ft. Finally, the well with the lowest correlation (well 5862603) has the most shallow screening interval at 100-152 ft. below the top of the casing and is the second closest well to the river. 43 Aquifer Heads and River Discharge 0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000 50,000 5/21/81 5/31/81 6/10/81 6/20/81 6/30/81 Date River DIscharge(cfs) 305 310 315 320 325 330 335 340 Aquifer Head (ft) River 5854801 5854706 5862603 Figure 4.5 Comparison of Daily River Discharge at USGS Gage 8159200 and Aquifer Heads at Three Nearby Wells This analysis reveals the difficulty in using sparsely recorded observational aquifer data to assess the interaction between the river and aquifer systems. Despite, the appearance of some correlation between the river and aquifer in observation well 5854706, with the lack of data taken at smaller time intervals, it was difficult to determine the response time of the aquifer to the changes in river discharge. Furthermore, the spatial variation in the aquifer?s response could not be assessed with the given data, because of the few number of wells available. There is also the matter of arbitrarily choosing a monthly average 44 of the daily streamflow values, which may have affected the correlation between the surface and groundwater fluctuations. As a result of this study, the direction of the research changed to include two linked physically distributed models of the river and aquifer systems. Unlike the observational data, a physically distributed model provides information in locations where groundwater data is not available. Furthermore, although a model may require additional field data, it inherently reduces the dependence on time series data such as the aquifer heads, which were not available. 4.3 Manning's Equation The physically distributed model of the Colorado River is based on Manning's equation and kinematic wave theory. In this study, Manning's equation, along with additional assumptions, was used to convert the river's flowrate, as determined by kinematic wave theory, into river stages. The stage was then added to the riverbed elevations to result in the river hydraulic heads, HR, for input into MODFLOW?s River Package. The Manning equation describes open channel flow based on three different parameters: hydraulic radius, R; friction slope, Sf; and n, the Manning roughness coefficient. The hydraulic radius R is defined as the cross-sectional area A of the channel divided by P, the wetted perimeter. The wetted perimeter is defined as the perimeter length of the channel that is in contact with the water, and thus "wetted." Sf is equivalent to the head loss along the channel, divided by the length of the channel. 45 In SI units, Manning's equation takes the following form: n SRV f 2/13/2= (4.3.1) where, V = velocity [L/T] R = hydraulic radius [L] Sf = friction slope n = Manning roughness coefficient In English units the above equation is multiplied by a factor of 1.49. SI units were used in constructing the surface water model. The Manning equation is only valid if the flow is turbulent. The criterion for turbulent flow is the following (Chow et al. 1988): 136 101.1 ??? fRSn with R in meters (4.3.2) The four additional assumptions used in this study are the following: ? Channel is a rectangular conduit 46 ? Width of the rectangular conduit is much larger than the depth of the water; wetted perimeter, P equals width of channel, w3 ? Manning roughness coefficient, n = 0.04 for clean, winding rivers (Chow et al. 1988) Solving Manning's equation for the river's water depth (y), with the aforesaid assumptions yields: 5/3 2/1 ? ? ? ? ??? ?= fwS Qny (4.3.3) where, Q = river flowrate (m3/s) w = width of the river (m) Sf = friction slope 4.4 Kinematic Wave Theory Kinematic wave theory describes the change in river flowrates, both in time and distance along the channel, based on lateral inflow into the river. In this study, the output of MODFLOW, namely the river leakage rates, is one component of the lateral inflow that is modeled. Therefore, the surface water 3 Future research may wish to eliminate this assumption by using the Newton-Raphson method to solve for R, the hydraulic radius. The Newton-Raphson method solves linear non-algebraic equations in an iterative fashion. 47 model uses kinematic wave theory and Manning?s equation in succession to determine the change in river hydraulic heads based on the aquifer lateral inflow. The kinematic wave model is a simplified version of the Saint-Venant equations that describe unsteady non-uniform flow in a channel. The entire set of the conservation form of the Saint-Venant equations is as follows (Chow et al. 1988): Continuity Equation qtAxQ =??+?? (4.4.1) Momentum Equation 0)(11 2 =????+?? ? ? ??? ? ? ?+ ? ? fo SSgx yg A Q xAt Q A (4.4.2) where, A = cross-sectional area of channel [L2] Q = flowrate of water in channel [L3/T] x = length along the river centerline [L] g = gravitational constant = [L/T2] t = time [T] y = water depth in river [L] q= lateral inflow into the river [L2/T] 48 The following are the assumptions of the Saint-Venant Equations1 as reported by Chow, et al. in 1988: ? The flow is one-dimensional; depth and velocity vary only in the longitudinal direction of the channel. This implies that the velocity is constant and the water surface is horizontal across any section perpendicular to the longitudinal axis. ? Flow is assumed to vary gradually along the channel so that the hydrostatic pressure prevails and vertical accelerations can be neglected. ? The longitudinal axis of the channel is approximated as a straight line. ? The bottom slope of the channel is small and the channel bed is fixed; that is, the effects of scour and deposition are negligible. ? Resistance coefficients for steady uniform turbulent flow are applicable so that relationships such as the Manning's equation can be used to describe resistance effects. ? The fluid is incompressible and of constant density throughout the flow. In the kinematic wave model, inertial and pressure forces are assumed to be negligible. With these assumptions, the kinematic wave equations are the following: 1 Neglecting the effects of wind shear, lateral inflow, and eddy losses and assuming the momentum correction factor, ? = 1 49 Continuity Equation qtAxQ =??+?? (4.4.3) Momentum Equation fo SS = (4.4.4) By setting the pressure and inertial terms of the momentum equation (Eq. 4.4.2) to zero and solving for A in terms of Q, Equation 4.4.4 could also be written as the following: A = ?Q? (4.4.5) where, ?,? = general coefficients Although this is a general solution, Manning's equation can be shown to be a specific solution of the momentum equation. Manning's equation (Eq. 4.3.1) along with the momentum equation (Eq. 4.4.4), yield the following relationships for a and b: 50 ??? ? ??? ?= oS nP 3/2a (4.4.6a) ? =0.6 (4.4.6b) Substituting Equation 4.4.5 for A into Equation 4.4.3 yields the kinematic wave equation: qtQQxQ =??+?? ?1bab (4.4.7) In this study there are two components of lateral inflow. These are qa and qr, representing the aquifer inflow and local surface inflow, respectively. Therefore, the form of Equation 4.4.7 that is used in the surface water model is shown in Equation 4.4.8. ra qqt QQ x Q += ? ?+ ? ? ?1bab (4.4.8) 4.5 Explicit Solution to the Kinematic Wave Equation Due to the complexity of the region being modeled, a finite difference solution was chosen over an analytical solution to the kinematic wave equation. The finite difference method is based on a Taylor series expansion of a differential equation, in this case Equation 4.4.8. The result of the finite difference method as applied to the kinematic wave equation is seen in Equation 51 4.5.1. This equation has been modified from what was reported by Chow et al. to incorporate both lateral inflow terms. Also note that the aquifer inflow term is held constant throughout the duration of the surface water run-time and therefore does not require time averaging. ? ? ? ? ? ? ? ? ??? ? ??? ? ++ ? ? ? ? ? ? ? ? ? ? ??? ? ??? ? ++?+ ??? ? ??? ? ++ ? ? = ? + + + + + + ?+ + + + + + 11 1 1 1 1 1 11 1 1 1 1 1 2 22 b b ab ba j i j i ai j ir j ir j i j ij i j i j i QQ x t qqqtQQQQxt Q (4.5.1) where, j = time-step index i = river reach index qr = overland runoff [L2/T] qa = aquifer recharge [L2/T] ?t = time-step length [T] ?x =river reach length (L) Two methods of solution for a finite difference equation exist: an implicit scheme and an explicit scheme. The implicit scheme involves solving the unknown values of the equation simultaneously and requires a specially constructed computer program. Alternatively, the explicit scheme solves the unknown values sequentially, in this case, with the help of an Excel spreadsheet. To ensure a stable solution, however, the explicit scheme requires smaller surface 52 water time-steps (Chow et al. 1988). In the explicit scheme, the length of the surface water time-step is constrained by the Courant condition described in Section 4.5.1. For the purpose of this study, the explicit solution was chosen for ease of implementation. 4.5.1 COURANT CONDITION A necessary but insufficient condition for stability of the explicit scheme is dictated by the Courant condition. This condition states that the time-step chosen for the surface water model must be smaller than or equal to the time it takes the wave to transverse any given river reach. Mathematically, the Courant condition is equivalent to Equation 4.5.2. k i c xt ??? (4.5.2) where, ?t = surface water model time-step [T] ?xi = reach length i [L] ck= kinematic wave celerity [L/T] The wave celerity, ck, can be determined from the Method of Characteristic solution of Equation 4.4.7. The Method of Characteristics produces a set of ordinary differential equations that are mathematically equivalent from a partial differential equation, in this case the kinematic wave equation. The Method of Characteristics solution to the kinematic wave equation 53 is presented in Equation 4.5.3. The right hand side equation can be rearranged to result in Equation 4.5.4, the wave celerity. By substituting in Manning's equation (Eq. 4.3.3), the wave celerity can be calculated using river parameters as shown in Equation 4.5.5. q dQ Q dtdx == ?1bab (4.5.3) 1 1 ?== babQcdt dx k (4.5.4) VnySc ok 3535 3/22/1 =?? ? ? ??? ?? ? ?? ? ?= (4.5.5) where, ck = wave celerity (m/s) Q = river discharge (m3/s) q = lateral inflow (m2/s) ?, ? = general coefficients So= riverbed slope y = water depth (m) n = Manning's roughness coefficient V = velocity (m/s) 54 4.6 Determining Initial Flowrate in the Colorado River In order to begin modeling the Colorado River, it was necessary to determine the initial condition of the river for the model. It was decided to use the average discharges in the Colorado as the initial condition. To determine the Colorado River?s average discharges, terrain and runoff data were modeled in a Geographic Information System (GIS). Because of the large amount of data involved in analyzing the entire Colorado watershed, a subwatershed in the domain of interest was the only one modeled in GIS. Contributions from the watersheds upstream of the model domain were determined by using average daily flowrate data collected from USGS stream gages. Section 4.6.2 discusses the GIS processing of the domain subwatershed. Section 4.6.3 describes the contribution of the upstream subwatersheds to the average flowrate. 4.6.1 DATA The data required to perform the average flowrate analysis is shown in Table 4.1. Although some of the data were in files at the Center for Research of Water Resources (CRWR), all files that were needed could just as easily have been downloaded from the Internet. All GIS processing of data for this study was conducted in version 3.2 of ArcView and version 7.1 of Arc/ Info. A Digital Elevation Model (DEM) is a digitized representation of the landscape and can be imported into ArcView as a grid. The value in each grid cell corresponds to the elevation of the terrain at the center of the grid cell. In this study, the grid cells were approximately 90 meters on each side. 55 The runoff grid was based on a coverage created by Andrew Romanek (Romanek 1998). The grid was created from data gathered by the USGS between 1951 and 1980. Table 4.1 GIS Data and their Sources Data Source Type of GIS Coverage 90 m. Digital Elevation Model (DEM) USGS Grid 8-Digit Hydrologic Unit Code Subwatersheds (HUC8) CRWR Polygon Counties CRWR Polygon River Reach 1 Files (Rf1) CRWR Polyline Runoff Grid CRWR Grid Width of Colorado River BEG -- Average daily flowrates upstream of domain watersheds USGS -- 4.6.2 GIS PROCESSING OF DATA To determine the contribution of the domain subwatersheds to the initial flowrate of the section of the Colorado River used in this study, it was necessary to cut the DEM into the shape of the domain subwatersheds. Subsequently, the DEM was taken through the steps of the CRWR-PrePro Project (Olivera 1998). From this process, a flow direction grid was produced. The runoff grid and the flow direction grid were then processed using the Non-Point-Source Project (nonpoint.apr) created in ArcView to produce a weighted flow accumulation grid. 56 Merging and Clipping DEMs To clip out the DEMs in the shape of the domain subwatersheds, the following steps were taken: ? Initially, Milam, Lee and Bastrop counties were ?intersected? in ArcView with the HUC subwatershed shapefile to determine which subwatersheds were underlying the three counties. Only one HUC8 subwatershed of the Colorado River, the Lower Colorado-Cummins (Cataloging Unit: 12090301), is located within these three counties. However, before the scope of this research was limited to modeling the Colorado River alone, the original GIS analysis was performed on all of the subwatersheds underlying the modeling domain. ? The 90 meter DEM data is available on the USGS web site (http://edcwww.cr.usgs.gov/doc/edchome/ndcdb/ndcdb.html). The DEM grid was then clipped out in the shape of the subwatershed basins of interest in Arc/Info. ? The DEMs were converted from their original floating point format to integer format to save disk space and decrease computing time in the subsequent analysis stage. Flow Direction Grid The DEM, now in the shape of the domain subwatersheds, was manipulated through the steps offered by CRWR-PrePro to create a flow direction 57 grid. First, the Rf1 line representing the Colorado and Brazos River were selected and clipped in the shape of the domain subwatersheds using the GeoProcessing Wizard Extension. Although more accurate digital representations of the Colorado and Brazos Rivers exist (in both River Reach 3 files and National Hydrologic Dataset formats), the Rf1 files are sufficient for the purpose of this study. The Rf1 representations of the Colorado and Brazos Rivers were then "burned into" the DEM by using the ?Burn-in Stream? command in the CRWR- PrePro menu. This command consists of raising all non-river grid cells by 5000 meters to ensure that water will flow into the location where the river is known to be. Due to inaccuracies in the DEM, depressions (i.e. grid cells bordered by two other grid cells with higher DEM values) can exist in the grid cells below the Rf1 representation of the rivers. To prevent water from stagnating in depressions in the river and disturbing the flow to downstream cells, the DEM was processed using CRWR-PrePro ?Fill DEM? function. During this process, depressions along the river were raised to the value of the lowest neighboring cell so as not to act like a sink in the river. Once the DEM had been filled, the flow direction grid was computed using the ?Flow Direction? command in the CRWR-PrePro Project in ArcView. In executing this command, the computer uses the 8-direction pour point model to calculate which direction the water will flow out of a cell. This model compares the DEM values of the 8 cells that surround any given cell and determines the direction of the steepest decline from the center cell. The model then assigns a number (and later a color) to this direction, based on the values shown in Figure 58 4.6. For example, if the direction of the steepest decline is due east, the number "1" is assigned to the center cell shown in the figure. In the resulting flow direction grid, seen in Figure 4.7, each color designates the direction of steepest descent and thus the direction that water will flow out of the cell. Although the flow direction grid contains data for all of the subwatersheds in the modeling domain, only the Lower Colorado-Cummins subwatershed was used to determine the initial flowrate in the Colorado River. 1 248 16 32 64 128 Figure 4.6 Schematic of the 8-Direction Pour Point Model 59 Figure 4.7 The Flow Direction Grid for the Subwatersheds Underlying Milam, Lee and Bastrop Counties Weighted Flow Accumulation Grid To determine the amount of runoff that gathered into the river cells, it was necessary to enter the flow and runoff grids into the ArcView Non- Point Source Analyst Project and employ the "Flow Accumulation" command. Unlike the flow accumulation command under CRWR-PrePro, which assumes each grid cell has an equal weight (i.e. 1), the Non-Point Source Analyst weighs the upstream cells with the runoff grid values before summing them up. For example, if 3 cells with runoff grid values of 4, 2 and 5 inches flow into a certain cell, the resulting flow accumulation value in the cell would be "11" inches and not "3" cells as reported by CRWR-PrePro. The flow eventually accumulates into the lowest DEM grid 60 cells, (i.e. the "burned in" river) and in this way calculates the average flow in the river being modeled. 4.6.3 UPSTREAM FLOWRATE Because the domain of interest was on the lower portion of the extensive Colorado watershed, it was necessary to consider the contribution of the flow originating from upstream of the domain subwatershed. Upstream flowrates were determined using the average daily flowrate provided by 3 USGS gages located directly upstream of the Lower Colorado - Cummins subwatershed. Figure 4.8 shows the location of the three gages with respect to the Lower Colorado - Cummins subwatershed. All three gages are located in Austin, TX. Gages 818600 and 815900 measure the flowrates of two tributaries of the Colorado River. They are located on Walnut Creek at Webberville Rd. and Onion Creek at U.S. Highway 183, respectively. Gage 815800 measures the river discharge of the main stem of the Colorado River at Austin. Table 4.2 shows the different dates of data used and the resulting flowrate at each gage station. Note that some days of data were missing from the interval stated in the third column of the table, accounting for the variation in the "Numbers of Days of Data" column from what is expected. The flowrate from the upstream watersheds of the Colorado River was determined to be 58 m3/s. 61 Table 4.2 Data from USGS Gages Located Upstream of the Subject Area USGS Station ID River/ Tributary Dates of Data Number of Days of Data Avg. Flowrate (m3/s) 8158600 Walnut Creek 12/29/74- 12/21/98 8758 0.9 8158000 Colorado River 12/29/74- 12/29/98 8694 54.4 8159000 Onion Creek 3/23/76- 10/17/98 9313 2.6 Figure 4.8 Location of the Three USGS Gages Used in Determining the Upstream Flowrate in the Colorado River 62 4.6.4 RUNOFF AND RIVER REACH LENGTHS In the surface water model, there are 35 reaches that make up the stretch of the Colorado River being modeled. These reaches correspond to the portions of the river overlying the River Package cells of the Colorado River in the 1st layer of the revised BEG MODFLOW model. As shown in Figure 4.9, the finite difference nodes defining the surface water model are located at the top and bottom of the river reaches, and at the boundary between the River Package cells, shown in the figure as blue boxes. Consequently, the most upstream MODFLOW River Package cell is located between nodes 0 and 1 of the surface water model. The river reach lengths, corresponding to Dx in Equation 4.5.1, were determined using the measuring tool in ArcView. The river flowrates resulting from the domain subwatershed were determined for each of the nodes 0-35, by querying the weighted flow accumulation grid at the node locations. This value was then added to the upstream flowrate and resulted in the initial flowrates for the Colorado River. Runoff values, corresponding to qr in Equation 4.4.8, were calculated based on the initial flowrate and river reach lengths according to Equation 4.6.1. Values of the river reach length and runoff can be viewed in Appendix A. ??? ? ??? ? ? ?= ? i ii ri x QQq 1 (4.6.1) 63 where, qri = overland inflow at node i [L2/T] Qi = initial flowrate at node i [L3/T] ?xi = length of river reach between nodes i and i-1 [L] Figure 4.9 Finite Difference Nodes Along the Colorado River 4.6.5 SURFACE WATER MODEL CALIBRATION The surface water model was calibrated based on the mean daily streamflow data of USGS gage 8160400 located on the Colorado River Above La Grange, TX. As shown in Figure 4.10, USGS gage 8160400 is situated near node 34. Streamflows from this gage were averaged over a span of 10 years starting in January 1, 1988. The 10-year average river discharge (81.65 m3/s) is approximately 9.4 m3/s more than what was predicted by the GIS analysis to be 64 the river discharge in node 34. To calibrate the model to match the observed data, the flowrate at node 0 was increased from 65.7 m3/s to 77.3 m3/s. New river heads were determined from these newly calibrated river discharge values, which caused a recalibration of the revised BEG model. To remain consistent with the original river heads in the BEG model, the riverbed elevations were adjusted to account for the new river stages, so that the revised BEG model contained approximately the same river heads as the original BEG model. Figure 4.10 Location of USGS Gage 8160400 with Respect to the Finite Difference Nodes 65 4.7 Excel-Based Surface Water Model To determine the effect of time-step on the interaction between a river and aquifer, a flood wave was simulated through the Colorado River. The rapid fluctuations in river discharge caused by the flood wave were averaged over various time-steps before being entered into the dynamically linked surface and ground water models. The organization of the surface water model was designed with the flood wave simulation in mind. Therefore, a description of the simulation is provided in Section 4.7.1 before the framework of the surface water model is discussed in the subsequent section. 4.7.1 SURFACE WATER SIMULATION The surface water model simulations were based on changing the inlet river discharge at node 0 to model the river discharges created by a flood wave. Each simulation occurs over a 32-day period. The inlet river discharge for the 32 days of the simulations were based on average daily streamflow data gathered by the USGS starting on January 8, 1991 from gage 8160400. This stream gage was also used to calibrate the surface water model and its location with respect to the domain can be seen in Figure 4.10. Although the gage is located near node 34, the streamflow data was applied as the boundary condition at the inlet of the domain (i.e. node 0). This was justifiable because the focus of this research is on the effect of time-step in modeling rapidly fluctuating river discharges; a baseline change of 4 m3/s in the river discharge is not expected to affect the results of this study. 66 Figure 4.11 presents the original USGS data for the inlet river discharge as 1-day averaged data. The inlet wave has three peaks occurring approximately during day 3, 11 and 28. The three peak discharges are approximately 1033 m3/s, 583 m3/s, and 427 m3/s. The river discharge in the long stretch of time between peaks 2 and 3 varies between 17 and 80 m3/s. To determine the effect of the time-step on modeling the river and aquifer interaction, 4 different time-steps were considered. The time-steps corresponded to averaging the river discharge over a period of 2, 4, 8 and 16 days. Figure 4.11 displays the time-step averaged river discharge at the inlet. As was expected, the smallest time-step was still able to capture the rapid fluctuations in the river discharge while the 16-day time step was not able to do so. 67 Inlet River Discharge 0 200 400 600 800 1000 0 8 16 24 32 Simulation Time (days) River Discharge (m 3 /s) 1-Day 2-Day 4-Day 8-Day 16-Day Figure 4.11 Time-Step Averaged River Discharge at the Inlet To ensure that any changes in the groundwater system were due to the inlet wave and not due to the new average river heads resulting from the surface water calibration, the groundwater and surface water systems were allowed to come into equilibrium. The dynamically linked models were run for 1280 days (or approximately 3.5 years) using 64-day time-steps. Although heads in the aquifer were still changing slowly with time, most of the changes occurred over a period that was longer than 32 days. The dynamically linked models were believed to be in pseudo-equilibrium for the purposes of this study. The new river and aquifer heads, now in pseudo-equilibrium, were used as the initial conditions for the 32-day simulations 68 There are three main components to the output of the simulations: river discharge, river hydraulic head and aquifer hydraulic heads. The output of the surface water model is comprised of the river discharge and hydraulic head in the 35 reaches. In order to assess the spatial variability in the aquifer's response to the flood wave, aquifer heads in two cross sections of the aquifer were recorded in the output. Figure 4.12 displays the location of the two cross sections in Layer 1. The first cross section includes nine cells in row 40, between columns 5 and 11. These cells include the aquifer cell directly below River Reach 14, shown as the light blue cell in Cross Section 1, and the four cells adjacent to it on both sides. The aquifer heads in Layer 2 of the same 9 cells were also recorded, so as to facilitate an understanding of the aquifer's vertical response to the flood wave. Cross Section 2 includes the ten aquifer cells located in row 37, between columns 8 and 17. These cells include River Reaches 27 and 28 along with the 8 cells in the same row straddling those reaches. Similarly, Cross Section 2 also extends to Layer 2. 69 Figure 4.12 Aquifer Cross Sections 1 and 2 in Layer 1 of the Revised BEG Model 4.7.2 SURFACE WATER MODEL The surface water model is an Excel file, ?kinematic8.xls?, which consists of 13 different worksheets. Four of the worksheets (inlet2, inlet4, inlet8, and inlet16) contain the time-averaged inlet river discharges for the four different time-steps. The kinematic wave and Manning's equations calculations are performed using six worksheets (Q, h, qa, qr, Par and qa,h_fix). The remaining three worksheets (flowrate, riv_head, and aq_head) are used to store the output of the dynamically linked models. The following sections provide a description of the methodology used in determining the surface water time-step as well as the worksheets in "kinematic8.xls." 70 Surface Water Time-Step To ensure that the Courant condition was met throughout the surface water simulation, a critical time-step was determined that was applied to the entire simulation. By examining the Courant condition equation (Equation 4.5.2) the critical time-step results from minimizing the value on the right hand side of the equation or the river reach length divided by the wave celerity. The critical time- step would need to satisfy the Courant condition in the smallest river reach and during the time of the greatest wave celerity. From Manning's equation, it is evident that wave celerity increases with increasing river discharge. Therefore, the first peak discharge of 1033 m3/s produces the highest wave celerity during the 32-day simulation. During that peak discharge, the time-step is also limited by the smallest reach, in this case Reach 15. A steady-state simulation in worksheet Q using 1033 m3/s as the inlet river discharge provided a critical time- step in Reach 15 of 197 seconds. A safety factor of 17 seconds was used in case the lateral inflow caused a significant increase in flowrate in any downstream reaches. Therefore, the critical surface water time-step used throughout each simulation is 180 seconds. Inlet Discharge Worksheets The time-step averaged inlet discharges are held in four different worksheets. The worksheets inlet2, inlet4, inlet8 and inlet16 contain the inlet river flowrates averaged over a time-step of 2, 4, 8 and 16 days, respectively. Although Figure 4.13 shows only a portion of inlet2, the other three inlet river discharge worksheets are organized in a similar fashion. 71 Starting from the far left, columns A-C contain different time variables used in the simulation. Column B of the worksheet contains the surface water time-step (180 seconds) and can be altered should other flood wave simulations be considered in the future. The inlet river discharge for each day of the simulation is located in each column, starting in column D. Therefore, in worksheet inlet2, columns D and E contain the 2-day average inlet flowrates for the first two days of the simulation. Similarly, the 2-day averaged inlet flowrates between days 2 and 4 are stored in columns F and G. The rows in the worksheet, starting at row 3, correspond to different surface water time-steps. There are 480 rows of flowrates in the worksheets, which correspond to the 480 surface water time-steps that comprise a day of the simulation. Although the inlet river discharge values did not change with each surface water time-step, the worksheets are set up to facilitate any kind of surface water data. Real-time data for instance, should it become available, could make better use of the 480 rows available in the inlet discharge worksheets. 72 Figure 4.13 A Section of the inlet2 Worksheet Kinematic Wave Worksheets The river parameters, such as riverbed slope, river width (or P, the wetted perimeter), and river segment length (?x), are contained in worksheet Par shown in Figure 4.14. RECNO corresponds to the record number given to the cell polygons in the GIS representation of the MODFLOW grid. Unlike, the inlet river discharge worksheets, the bold numbers along the top of the 4th row correspond to each of the 35 river nodes. The most important parameter on this worksheet is the lumped parameter ?, one of the main parameters in Equation 4.4.8. The Manning coefficient factor used throughout the simulation is 0.04, for clean winding rivers, and can be seen in cell B2. In addition, the Par worksheet contains the river leakage rates and the river discharge in the 35 reaches that resulted from the calibration of the dynamically linked models. These values are 73 used to produce consistent initial conditions in the surface water model for each simulation. Figure 4.14 A Section of the Par Worksheet Like the Par worksheet, qa and qr provide input values for the surface water model. Worksheets qa and qr contain values of aquifer recharge and runoff, respectively. In the qa worksheet, because the river recharge by the aquifer is assumed to stay constant for the entire surface water simulation, there is only one row of aquifer recharge values. The qa values are based on the output of MODFLOW, but have been modified in the qa.h_fix worksheet to stay consistent with the units used in the surface water model. The qr worksheet, on the other hand, contains a row of runoff data per surface water time-step. It was believed at the beginning of the project that precipitation events might also be considered when modeling the surface water system. These events would give rise to time- dependent runoff, and the worksheet was set up to accommodate this aspect of the 74 analysis. However, due to the added complexity posed by simulating precipitation, simulations that were conducted in this study were limited to the steady-state average runoff values that were calculated in Section 4.6.4. The qr worksheet, however, is equipped to handle the analysis of time-dependent runoff, should further analysis of the region be desired. The qa,h_fix worksheet exists to modify and reorder the two parameters, HR and qa, that are being transferred between the two models. The parameters in the surface water model are ordered from upstream to downstream. The revised BEG MODFLOW model, however, is organized by increasing row number, which at times does not correspond to the direction of river flow. The sequence of the qa and HR vectors are therefore shuffled in qa,h_fix to suit the organization of the model that is being updated. In addition, qa,h_fix worksheet converts MODFLOW?s output (i.e. river leakage rates in units of ft3/d) to suit the units of qa (m2/s) in the surface water model. Lastly, in qa,h_fix the riverbed elevations are added to the river stage values to result in river hydraulic head before MODFLOW?s River Package was updated. Worksheets Q and h are set up very similarly to each other and contain what can be considered the output of the surface water model. The Q worksheet contains the river discharge at the 36 finite difference nodes, which are calculated using the linear explicit scheme described in Section 4.5. Similar to the Par worksheet, the columns of worksheet Q contain the river discharge by each river node shown in row 4. The inlet (or node 0) river discharge is stored in column D. During each day of the simulation, depending on the time-step specified by the 75 user, a column from one of the four inlet river discharges worksheets is copied into column D of the Q worksheet to provide the new boundary condition for that day of the simulation. Changes in the new inlet discharge create a wave that is propagated through the rest of the reaches via the kinematic wave equation. The h worksheet is directly linked to the Q and Par worksheets and calculates the river stages based on Manning?s equation. Figure 4.15 A Section of the Q Worksheet Output Worksheets Three worksheets contain the output of the dynamically linked surface and groundwater models. The river discharge and river hydraulic head are recorded in the riv_head and flowrate worksheets, respectively. The aq_heads worksheet contains the aquifer heads from in Cross Sections 1 and 2, as well as the river leakage rates for Reaches 14, 27 and 28 and the total river leakage for all 35 reaches being modeled. 76 CHAPTER 5: VISUAL BASIC INTERFACE 5.1 Introduction The Visual Basic interface (called the ?Interface? for the remainder of this document) was created inside the Visual Basic editor of Excel. The purpose of the Interface is to act as the dynamic link between the MODFLOW model and the Excel-based surface water model. More specifically, the Interface elicits input from the user, transfers data between the two models and executes each model in a cyclic fashion. The Interface also writes the river hydraulic heads, the aquifer hydraulic heads in the two cross sections of the aquifer, and the river discharge to three output worksheets in the ?Kinematic8.xls? Excel file. Visual Basic was chosen as the programming language of the interface because of its compatibility with Excel; Excel uses Visual Basic to perform many of its commands. Section 5.2 gives an overview of the capabilities of the Interface. In Section 5.3, the basic structure of the Interface?s programming code is described. Section 5.4 discusses some of the format requirements for the MODFLOW packages. In conclusion, possible methods to improve the interface performance are discussed in Section 5.5. 77 5.2 Overview On a very general level, the function of the Interface can be broken down into three main objectives, as depicted in Figure 5.1: reading input, dynamically linking the models and writing output. Figure 5.1 Schematic of the Three Main Functions of the Interface 5.2.1 INPUT The input to the Interface is entered by way of a Visual Basic form called ?frmInput?, shown in Figure 5.2. The variable names, shown in the Input box of Figure 5.1, store the user-defined values that are entered through "frmInput". In order to stay within programming conventions, single letters (such as the "d" in dTotTime) were used as prefixes to variable names to indicate the variable type 78 of each variable. The letter "d" and "s" designate double precision and string variable types, respectively. The first of the three user-specified variables is dTotTime, the total simulation time of the dynamically linked models. Although, river discharge data exists for 32 days of simulations, the length of the simulation is still kept as a user-defined variable for added flexibility. Should any simulations that involve more or less days of data be investigated, the user would be able to change this variable without having to change the Interface code. This feature was also extremely helpful during the debugging of the Interface code, when a number of shorter simulations were run to determine if certain portions of the model were working correctly. Next, ?frmInput? asks for the user to designate a time-step for the simulation via the variable dTimestep. Since the surface water model only contains values for four time-step averaged inlet discharge worksheets (inlet2, inlet4, inlet8 and inlet 16), the Interface only simulates flood waves for dTimestep equivalent to 2, 4, 8, and 16 days. All other time-steps result in groundwater/surface water simulations under non-flood wave conditions. Finally, the input form asks the user for the pathname of a folder ?modflow? that contains the following: the surface water model (?Kinematic8.xls?), the MODFLOW-96 executable and the individual MODFLOW packages. The pathname of the folder entered by the user is stored in the variable sPath. 79 Figure 5.2 "frmInput" Asks the User for Input into the Visual Basic Interface 5.2.2 DYNAMIC LINKAGE The second capability of the interface is its main function, namely, dynamically linking the two models. This subject will be covered in much more detail in the following section, but will be simplified to provide a quick overview here. The dynamic linkage can be thought of as one large loop. The surface water model is run first with its updated aquifer inflow rates and inlet river discharge, followed by MODFLOW with its updated river hydraulic heads. After each time the MODFLOW is run, the output is written to the three output worksheets. The loop is repeated until the total simulation time entered in by the user, or dTotTime is reached. At this time the interface ends the Visual Basic program, leaving the output of the models stored in the three output worksheets. 80 The surface water model isn?t run in the conventional sense, but because the cells in the Q worksheet are linked to the cells in the aquifer lateral inflow (qa) worksheet, the updated qa and inlet discharge values automatically result in new river flowrates and new hydraulic heads. The MODFLOW program was run conventionally, although attempts to automate this function in the interface were unsuccessful due to compiler incompatibility. This topic is discussed in more detail in Section 5.5. 5.2.3 OUTPUT The last capability of the interface is the writing of the output of the MODFLOW and surface water models to the three output worksheets. The output worksheets are updated after each model is run. The river's discharge and hydraulic heads are fairly easily copied from the Q and qa,h_fix worksheets, and pasted into the flowrate and riv_head sheets, respectively. The MODFLOW output, however, is more difficult to extract. Manipulating MODFLOW?s listing file, where all printed output is written, has two challenges. As described in Chapter 3, depending on what was indicated in the Output Control package, the length of MODFLOW's text output or listing file can vary. For instance, the Output Control Package can determine whether the starting aquifer heads are printed or not printed to the listing file, thus changing the location of the aquifer heads that are to be extracted within the file. To allow for flexibility in the listing file format, the ?Find? capabilities of Excel were employed to locate key words in the listing file, which would indicate the location of the aquifer heads of interest. 81 The aquifer heads cells were then referenced from the row containing the key words. The second challenge in extracting the MODFLOW aquifer heads of the Colorado River Package cells, lay in the format of the output file, known as the listing file. MODFLOW -96 formats the listing file so that each row in the domain would fit onto an 8 1/2? by 11? paper when printed. Therefore, though River Package cells were next to each other in the domain, they were rows apart in the listing file. To overcome this problem, a sub-procedure named ?makegrid? was written to convert the output file back into the domain grid format. This sub- procedure, along with others, is described in the following section. 5.3 Structure The Interface contains two main procedures and nine sub-procedures. The two main procedures result from clicking either the OK or Exit buttons located on bottom of the ?frmInput? form. If the Exit button is clicked, the ?cmdExit_click? procedure is called, which ends the Visual Basic program. Alternatively, clicking the OK button calls the ?cmdOK_click? procedure that proceeds to run the dynamic link. To modularize the dynamic linking process, nine sub-procedures were created to perform the separate tasks within the ?cmdOK_click? procedure. These sub- procedures are described in Table 5.1. Figure 5.3 provides a flow- chart to show how the sub-procedures are organized to perform the dynamic linkage. 82 Table 5.1 The Nine Sub-Procedures That Constitute the "cmdOK_Click" Procedure Sub-procedure Function Initialize Clears output sheets and initializes the river leakage rates and river discharges Surface_Run Sets new inlet river discharges as the upper boundary condition for one day simulation Outputsheet 1. Runs Surface_Run for the number of days in iTimestep 2. Copies the river flowrate and hydraulic heads into the flowrate and riv_head sheets Rename Copies the binary file containing the heads at the end of one MODFLOW run to the binary file containing the initial heads for the next MODFLOW run Open_Bas Enters the new time-step of the MODFLOW run into Basic Package Open_Riv Updates river hydraulic heads in River Package Run_MODFLOW Runs MODFLOW-96 executable Output_qa 1. Extracts aquifer heads from MODFLOW's listing file and pastes into the aq_head worksheet using "Makegrid" and "Outputhead" sub-procedures 2. Extracts river leakage rates for each of the 35 River package cells from MODFLOW's listing file Makegrid Reformats aquifer heads in the MODFLOW listing file to resemble the domain grid Outputhead Extracts aquifer heads underlying Colorado River from the MODFLOW listing file 83 The Visual Basic code corresponding to each sub-procedure can be viewed in Appendix B. Initialize Surface_Run Outputsheet Rename Open_Bas Open_Riv MODFLOW_Run Output_qa iLoopcount = iLoopcount + 1 Time = Time + iTimestep Time >= TotTime MODFLOW Run No End Yes Surface Water Model Run Figure 5.3 Flowchart of the "cmdOK_click" Procedure 84 5.4 Formatting of MODFLOW Files The formatting of the MODFLOW input files is crucial to the performance of the interface. Even with the free format option offered in the Basic Package, MODFLOW-96 can only read space and comma delimited files. The two files that are opened and saved in Excel during the coupled-model run are the Basic and River packages. Due to its large array of numbers, the Basic Package would not save properly with the space-delimited format and therefore was saved using the comma-delimited format. To stay consistent with the Basic Package formatting, the River Package was also saved as a comma-delimited file. The format of the MODFLOW files must be strictly adhered to because the Interface opens the file in Excel with these two formats incorporated in the code. 5.5 Interface Enhancement There are several changes that would enhance the performance of the Interface. These changes are described in the following sections. 5.5.1 AUTOMATION OF MODFLOW RUN The interface was unable to automate the process of running the MODFLOW executable. This limitation is in part due to the structure of MODFLOW's source code. Once MODFLOW is run by the Interface, a window emerges asking the user to enter in the name of the Name file of the model of interest. As described in Chapter 3, the Name file contains a list of the names of 85 the different packages that are used during the simulation. Attempts to automate the process of typing in the name of the revised BEG model's Name file (cw409_s2.nam) were unsuccessful. Therefore, the user is required to manually enter in the name of the Name file into MODFLOW after each time MODFLOW is run. One way of overcoming this problem would be to change MODFLOW?s source code to automatically search for the Name file and recompile. This approach was not taken because of the compiler-related errors that were anticipated to occur. The format of the files had been adjusted to be compatible with the Lahey compiler used by the USGS in its downloadable executable. However, the only FORTRAN compiler available for this work was the Microsoft Power Station FORTRAN compiler, which most likely would have been incompatible with some of the input files and function calls. 5.5.2 EFFICIENCY The surface water model portion of the dynamic linkage was significantly longer in execution time than the MODFLOW model. Two processes contributed to the relatively long computation time of the coupled programs. The first was the Excel-based explicit solution. Due to the Courant condition, the explicit solution requires smaller surface water time-steps than what is required by the implicit solution. These smaller time-steps, in turn, result in a longer surface water model run-time. The ?makegrid? sub-procedure also hindered the efficiency of the Interface. This procedure changes the format of MODFLOW?s output file to a 86 format resembling the domain grid. Although the procedure made it easier to extract the aquifer heads of interest and to subsequently check if the correct aquifer heads had been extracted, it significantly decreased the efficiency of the model. 87 CHAPTER 6: RESULTS 6.1 Introduction As described in Chapter 4, in order to investigate the aquifer?s response to changing river heads, a flood wave was simulated to pass through the modeling domain. The influence of time-step on the river and aquifer interaction was explored by averaging the river discharge at the upper boundary of the modeling domain for 2, 4, 8 and 16-day time-steps. The daily streamflow data from USGS gage 8160400 was averaged over the time-step prior to being used as input into the surface water model. The river heads produced at the end of the surface water model run were subsequently used to update MODFLOW?s River Package. The changes in river discharge and aquifer hydraulic heads, as a function of the time-step, are discussed in this chapter. Section 6.2 investigates the spatial variation in flowrate along the length of the river. Next, discussions of the changes in aquifer head in Layer 1 and Layer 2 are presented in Sections 6.3 and 6.4, respectively. Section 6.5 presents the trend observed in river leakage rates as a function of time-step. Finally, Section 6.6 discusses the appropriate approximation of the Saint-Venant equation, should flood waves not be considered in future simulations. 88 6.2 Spatial Variation in River Reaches As seen in Figure 6.1, the inlet discharge was the main component in dictating the flowrates throughout the 35 reaches. During each of the three 2-day time-steps shown in the figure, all 36 finite difference nodes experienced approximately the same river discharge. This is not surprising, considering the time it takes any change in inlet flowrate to transverse all 35 reaches ranges between approximately 0.5 and 2.0 days. These values were calculated by summing up the wave travel times for each river reach under the maximum and minimum inlet river discharges of the 2-day time-step simulations, respectively. Therefore, even the 2-day time-steps were often too long to capture any transitions the downstream reaches were experiencing due to the changes in inlet flowrate. 89 River Discharge Along River 110 120 130 140 150 160 170 0 20 40 60 80 Distance from Inlet (km) River Discharge (m 3/ s) Days 6-8 Days 10-12 Days 14-16 Figure 6.1 River Discharge Along the Colorado River During Three 2-Day Time-Steps The small increase in the discharge along the river is caused by the lateral inflow terms, qa and qr. The aquifer lateral inflow (qa) can be either positive or negative depending on the river flowrate, while the runoff lateral inflow (qr) is constant throughout the simulation and remain positive. A comparison of the summation of the qa and qr terms for the 35 reaches revealed that the range of the lateral aquifer inflow is between two to three magnitudes smaller than the runoff inflow. The runoff flow (qr) of the 74.4 km. stretch of the Colorado River being modeled was held at a constant rate of 8.3 E-5 m2/s for the total simulation time. The absolute average of the aquifer lateral inflow (qa) for that same stretch of river was 1.3 E-6 m2/s and 2.0 E-7 for the 2-day and 16-day time-steps, respectively. The average aquifer lateral inflow (qa) resulting from the steady- 90 state analysis conducted prior to the 32-day simulations was calculated to be 9.3 E-8 m2/s. The river head resulting from the qa induced changes in the flowrate are therefore significantly smaller when compared to those induced by the qr. Furthermore, both lateral inflow terms are significantly less influential than the changing inlet river discharge. Therefore, apart from providing insight to the amount of water that was transferred between the two systems, the aquifer lateral inflow had little effect on the river discharge. The overlapping graphs presented in Figure 6.2 show a different interpretation of the same results. The different river nodes in Cross Sections 1 and 2 corresponding to Reaches 14 and 27, shown in Figure 4.7, experienced almost identical changes in river discharge as a function of time. Consequently, similar trends were observed in the river heads and the aquifer heads in both cross sections. For the sake of brevity, only the results from the analysis of Cross Section 1 will be discussed in this chapter. However, it can be inferred that the all of the trends that pertain to Cross Section 1 were also observed in Cross Section 2. The results of the simulations for Cross Sections 1 and 2, as well as the Colorado River, are presented in Appendix C 91 Flood Wave Along the River 0 100 200 300 400 500 600 700 800 900 0 8 16 24 32 Simulation Time (Days) River Discharge (m 3 /s) Inlet Node 14 Node 27 Figure 6.2 River Discharge as Function of Simulation Time at Three Locations Along the River. 6.3 Aquifer Head in Layer 1 Although the simulations provided results for the eighteen cells comprising Cross Section 1, the results from six cells, in particular, are presented in this section. Figure 6.3 shows the location of the three cells in Layer 1 of Cross Section 1. The other three cells are located in Layer 2 of the MODFLOW model, as shown in Figure 6.4. All three cells have roughly the same width, ranging from 1 mile for the cell containing Reach 14 and its neighboring cell, to 1.25 miles for the cell two away from Reach 14. 92 Figure 6.3 Location of Selected Cells in Cross Section 1 Figure 6.4 3-Dimensional View of the Selected Cells in Each Layer of Cross Section 1 93 6.3.1 EFFECT OF TIME-STEP The amplitude of the fluctuations in the river and aquifer head is a function of the time-step used in the simulation. Figure 6.5 presents the changes in the aquifer head in the Reach 14 River Package cell for three of the four time-steps. During the 2-day time-step simulation, the aquifer head below Reach 14 increased to 108 ft. during the first peak discharge in the river. However, the results of the 16-day time-step simulation show that the aquifer in the same cell increased to only 88 ft. during that period. Therefore, fluctuations in both the river and aquifer heads were observed to decrease with an increasing time-step. This result can be attributed to the fact that larger time-steps average over the peak discharges and consequently, by Manning?s equation, the river hydraulic heads. Since the river head is the driving force for any changes in the aquifer in the vicinity of the river package cells, reductions in the variability of the river heads as a result of longer time-steps produce smaller head changes in the aquifer. 94 Aquifer Head Below Reach 14 70 75 80 85 90 95 100 105 110 0 8 16 24 32 Simulation Time (days) Aquifer Head (ft) 2-Day 4-Day 16-Day Figure 6.5 Changes in Aquifer Head in the Cell Containing Reach 14 6.3.2 EFFECT OF DISTANCE FROM RIVER The spatial variability in the aquifer's response to the flood wave was examined by analyzing the aquifer heads in Cross Section 1. Figure 6.5 presents the aquifer heads corresponding to the three cells in Layer 1. Although, not shown here, the changes in river head in Reach 14 as calculated by the surface water model, is identical to the aquifer head below Reach 14. Therefore, during each interface time-step, the MODFLOW model allowed the aquifer head in the River Package cell to come into equilibrium with the newly calculated river head. It should be noted that this result is influenced by the hydraulic conductivity of the riverbed. The ease with which water is exchanged between the aquifer and river systems depends on the assumptions made regarding the texture of the 95 riverbed material. As noted in Chapter 3, the hydraulic conductivity chosen for this study was determined to be characteristic of silt or loess soils and therefore deemed suitable. However, since assumptions were made regarding the hydraulic conductivity and the riverbed thickness, both of which contribute to the riverbed conductance, a sensitivity analysis of the river conductance term is recommended for any future research of this system. Heads in Layer 1 of Cross Section 1 70 75 80 85 90 95 100 105 110 0 8 16 24 32 Simulation Time (days) Aquifer Head (ft) Reach 14 One cell away Two cells away Figure 6.6 Aquifer Heads in Layer 1 of Cross Section 1 Using a 2-Day Time-Step Two observations can be made regarding the aquifer heads in Cross Section 1 as shown in Figure 6.6. The first is that the fluctuations in aquifer head diminish as the distance from the river increases. The aquifer cells neighboring Reach 14 show less dramatic changes in aquifer head as a result of the flood wave 96 in the river than those observed in the Reach 14 River Package cell. For the 2-day simulation, the aquifer head fluctuations in the Reach 14 cell and those in one cell and two cells away from Reach 14, show a range of approximately 30 ft., 7 ft., and 1 ft., respectively. Therefore, the wave in Layer 1 of the aquifer appears to dampen or attenuate as it flows away from the river. The second observation is that cells farther from the river display a lag in response time to the changes occurring in the river. Unlike the immediate response seen in the aquifer head in the cell containing the river reach, the head fluctuations in the neighboring cells were slower to respond to the oscillations in the river. Figure 6.6 shows that there was only one peak in the cell located two cells away from Reach 14, which occurred approximately 8 days after the initial peak in river discharge. The simulation time appears too short to fully capture the cell's response to the two subsequent peak discharges. In conclusion, aquifer heads situated farther from the river appear to fluctuate less and respond more slowly than those directly below the river. 6.4 Aquifer Head in Layer 2 The aquifer heads in Layer 2 shows a different response pattern than those in Layer 1. Changes in aquifer head in Layer 2 in the three selected cells of Cross Section 1 are presented in Figure 6.7. Due to the barrier imposed by the semi- impermeable layer separating Layers 1 and 2, the range in head fluctuation is considerably smaller than those observed in Layer 1. As shown in the figure, the head changes in Layer 2 are on the order of inches. In addition, unlike the heads 97 in Layer 1, the Layer 2 heads neighboring Reach 14 do not show any lag time in responding to the flood wave. It should be noted that there is a hint of a lag time seen in Layer 2, in the cells three and four cells away from the river reach. However, this lag time is much less dramatic than what is observed in Layer 1. Furthermore, the heads in Layer 2 display very little wave dampening away from the river. All three cells in Figure 6.7 show a similar range in head fluctuation. Heads in Layer 2 of Cross Section 1 101.5 101.6 101.7 101.8 101.9 102 102.1 102.2 102.3 102.4 102.5 0 8 16 24 32 Simulation Time (days) Aquifer Head (ft) Reach 14 One cell away Two cells away Figure 6.7 Aquifer Heads in Layer 2 of Cross Section 1 Using a 2-Day Time-Step 98 6.5 River Leakage Rates The river leakage rates for Reach 14 and the segment of the Colorado River in the modeling domain are presented in Figures 6.8 and 6.9. As seen in the figures, there is a positive correlation between river discharge and river leakage rate. The peak discharges correspond to the peak river leakage rates. This is evident from Darcy's Law in which higher river hydraulic heads produce higher leakage rates. The peak discharge values were significant enough to change the segment of the river being modeled from a gaining stream to a losing stream. The figures are presented in units of m3/s to provide a comparison of the leakage rates and the river discharge. Therefore, for a 1033 m3/s river discharge, the leakage rate for the 74.4 km stretch of the Colorado River was less than 0.50 m3/s. These results are based on the initial head in the aquifer, the hydraulic conductivity of both the riverbed and Layer 1, as well as other factors. 99 River Leakage Rate for Reach 14 -1.0E-03 -5.0E-04 0.0E+00 5.0E-04 1.0E-03 1.5E-03 0 8 16 24 32 Simulation Time (days) River Leakage Rate (m 3 /s) 2-Day 4-Day 16-Day Figure 6.8 River Leakage Rate in Reach 14 Using Three Different Time-Steps River Leakage Rate for Colorado River -0.300 -0.200 -0.100 0.000 0.100 0.200 0.300 0.400 0.500 0 8 16 24 32 Simulation Time (days) River Leakage Rate (m 3 /s) 2-Day 4-Day 16-Day Figure 6.9 Sum of River Leakage Rates for All 35 Reaches Using Three Different Time-Steps 100 Interestingly, the minimum river leakage rate did not result from the lowest river discharge. In fact, the 2-day and 4-day time-step simulations, show the lowest river leakage rates (i.e. the highest amount of baseflow to the river) occurred right after the peak discharges. Figure 6.10 shows the aquifer head change across Cross Section 1 in Layer 1 as a result of the flood wave. The x-axis of the figure is the distance in miles from Reach 14. Negative distances denote that the cells are located west of the river reach (or left in Figure 6.3). The cell containing Reach 14 is therefore at the intersection of the two axes. The head gradient under the initial condition shows that groundwater is typically recharging the Colorado River. During the first peak discharge (during days 2-4) the head gradient is inverted, causing the Reach 14 to change from a gaining stream to a losing stream. However, once the first wave had passed, (during days 4-6) the piezometric head profile resumed its original gradient towards the river. The negative river leakage rates following the peak discharge indicate that the water that was leaked from the river during the peak discharge returned to the river as baseflow due to the return of the piezometric head surface to its original profile. 101 Aquifer Heads in Layer 1 of Cross Section 1 20 40 60 80 100 120 -5 -4 -3 -2 -1 0 1 2 3 4 Lateral Distance from Reach 14 (miles) Aquifer Head (ft) Initial Days 4-6 Days 6-8 Figure 6.10 Aquifer Heads in Layer 1 of Cross Section 1 The total volume of surface water lost to the aquifer as leakage during the 32-day simulations is presented in Table 6.1. These values were calculated by summing the river leakage rates, as determined by MODFLOW, over the 32-day total simulation time. The amount of water exchanged was examined for three segments of the Colorado River: Reach 14, Reach 27 and the entire segment of Colorado River being modeled in this study. The table presents values for the volume of water exchanged between the two systems based on 2, 4, 8 and 16-day time-steps. The volume exchanged is reported in the same units used by the revised BEG MODFLOW model, namely cubic feet. Initially, it appeared as if for each simulation and averaging method, the minimum gain to the river occurred using a time-step of 4 days. This trend was 102 difficult to decipher. On further investigation, it became apparent that the third peak occurred right before the end of the 32-day simulation. Therefore, for time- steps greater than 2 days, it was possible that the response of the aquifer to the third peak was not able to be determined because the simulation came to an end. To rectify the possibility of the end of the simulation affecting the results, the sum of the river leakage rates for the 35 reaches was examined over the first 8, 16, 24 days of the simulations, before the occurrence of the third flood wave peak. These values are presented in Table 6.2 along with the results of the 32-day simulation. 103 Table 6.1 Volume of Groundwater Flow Entering into the Colorado River over the 32-Day Total Simulation Time Volume of River Leakage (ft3) Time-Step (days) Reach 14 Reach 27 Total 2 -10115 -3617 -691848 4 -7770 1011 -181978 8 -8016 64.5 -404740 16 -8368 -905 -576911 A comparison of the different values in Table 6.2 indicates that a pattern may exist between time-step and the volume of water exchanged between the river and aquifer. Although seemingly anti-intuitive, the results show that baseflow to the river increases (or river leakage rate decreases) with a decreasing time-step. Prior to obtaining these results, it had been hypothesized that the smaller time-steps would capture higher peak flows and therefore would result in an increase in river leakage. Assuming that the results accurately model what occurs in the aquifer, these results point toward the importance of the initial piezometric heads of the aquifer in affecting the river leakage. As described above, there is a large influx of groundwater into the river in the time-step after the peak discharge. The high baseflows that occur after the peak discharges, 104 which are the most evident during the 2-day time-step simulation, may be the factor that causes the trend in the volume of water leaked. Therefore, along with other parameters, the initial piezometric head profile may be worth investigating. If the initial head gradient was away from the river, possibly due to well pumping, water leaked from the river under high river discharge will most likely not return to the river. Therefore, the trend observed in the volume of water exchanged between the two systems under a losing stream scenario, may be different than what was observed in the gaining stream scenario presented here. Table 6.2 Volume of Water Leaked from River During the First 8, 16, 24 and 32 Days of Simulation Volume of Leakage from 35 Reaches (ft3) Simulation Time (days) Time-Step (days) 8 16 24 32 2 45652 -351852 -813655 -691848 4 119198 -48947 -765617 -181978 8 713439 328280 -657081 -404740 16 -- 461068 -- -678925 However, due to the preliminary nature of these results, it is important to determine whether what is being modeled is actually what is occurring in aquifers under similar stresses. One possible explanation of the results may be that there is either an error or inaccuracy in the modeling systems. This could be an error in 105 how the two models are linked, in which case the output of another coupled modeling system (such as MODFLOW-SURFACT or MIKE-SHE) would need to be provided for comparison. Another possibility is that the models are working properly, however, basic assumptions in either the MODFLOW or surface water model create an inaccuracy. This inaccuracy may be compounded by the number of model runs or may be exacerbated by high river discharge. Two possible assumptions that could be the cause of the inaccuracies are the use of Darcy's Law in modeling the river and aquifer interaction and the riverbed slopes used in Manning's Equation. Both of these assumptions fall in the latter category of becoming increasingly inaccurate with higher surface water flowrates. To determine if these assumptions are creating inaccurate results, these results would need to be compared to observational data gathered from wells near the river during a similar flood wave. Therefore, before the results are assumed correct, the dynamically linked models must be validated. 6.6 Approximation of the Saint-Venant Equations Future investigations of the aquifer river interaction in this study area may wish to include simulations that are not influenced by flood waves. In Chapter 2, the criteria presented by Vieira (1983) in determining the appropriate approximation to the Saint-Venant equations when modeling lateral inflow were discussed. In order to determine the correct wave routing method under a non- flood wave condition, the two governing parameters of the criteria, k and Fo were 106 calculated for the dynamically linked models used in this study. The equations pertaining to the two parameters are presented in Equations 2.3.1 and 2.3.2 and again here as 6.6.1 and 6.6.2. Equation 6.6.3 defines the Chezy coefficient in terms of Manning's roughness coefficient. 3/1 24 3 sin ??? ? ??? ?= qC Lgk q (6.6.1) 2/1tan ??? ? ??? ?= gCFo q (6.6.2) n RC 6/1= (6.6.3) where, k = kinematic wave number g = gravitational acceleration [L2/T] L = length of river reach [L] ? = constant angle of the slope R = hydraulic radius [L] q = lateral runoff [L/T] n = Manning's equation Fo = Froude number The initial conditions of the river, resulting from the coupled model calibration, were used in calculating the two governing parameters. Furthermore, 107 the 35 separate reaches were modeled as one unit. Therefore, the average initial stage in the river was used in calculating the hydraulic radius, R. The slope was determined by calculating the elevation drop between the upstream and downstream ends of the river, and by dividing by the river length. Due to the order of magnitude difference between the aquifer lateral inflow (qa) and the runoff lateral inflow terms (qr), qa was neglected when calculating both k and Fo. The resulting Fo and k values were 0.090 and 22.9, respectively. The flow in the river is therefore subcritical flow. According to the criteria published by Vieira, these values indicate that lateral inflow into a river under non-flood conditions, should be modeled using the diffusive approximation. Therefore, the diffusive wave equation will be required for any further studies of the aquifer and river interaction in the study area under a steady-state upper boundary condition. 108 CHAPTER 7: CONCLUSIONS 7.1 Introduction This chapter presents a summary of the results of the dynamic linkage of a groundwater and surface water model. Section 7.2 discusses the methodology used to meet the research objectives outlined in Chapter 1. Section 7.3 presents a discussion on the different software issues that became evident in linking the two models. A summary of the results of the time-step simulations is presented in Section 7.4, while recommendations for future research are presented in Section 7.5. 7.2 Objectives The following methodologies were used to meet the three research objectives discussed in Chapter 1: ? The first objective was to create a physically distributed surface water model of the Colorado River. The model consisted of 13 Excel worksheets and used the kinematic wave equation and Manning's equation in succession to determine the change in hydraulic river head caused by lateral inflow and a transient upper boundary condition. The worksheets in the model include worksheets containing input inlet river discharges, 109 river parameters, kinematic wave and Manning?s equation calculations as well as worksheets to store the output of the dynamically linked models. ? The second objective was to design the code for the interface that dynamically linked the models. In this project, the interface was written in the Visual Basic Language in the Excel Visual Basic editor. The interface was able to elicit input from the user regarding the time-step, run the surface water and groundwater models in an alternating fashion as well as write the output to worksheets in the Excel surface water model. Due to a lack of a Lahey FORTRAN compiler, the running of the MODFLOW model was not fully automated in this research. ? The effect of time-step in modeling the interaction between the groundwater and surface water systems was analyzed by simulating a flood wave through the Colorado River. The results of the simulation, however, are preliminary and need to be verified by a comparison with results produced by another surface water/ groundwater modeling system and with observed data in stream gages and aquifer wells located near the Colorado River. The preliminary results of the simulations are presented in Section 7.5. 110 7.3 Software In creating the dynamic linkage between the two models, there were observations made regarding the role and potential role for software in facilitating the dynamic linkage. Obstacles caused by using incompatible software were also encountered. These observations are discussed in the following subsections. 7.3.1 THE ROLE OF GIS A Geographic Information System (GIS) was used to assemble data for input into the surface water model, and can be useful in providing the data for MODFLOW's input files. Equally as important, the GIS provided a geographically referenced framework in which the surface and groundwater modeling domains could overlap. In this way, river reaches were created specifically to overlie MODFLOW's River Package cells. A GIS, therefore, allowed the interaction between the two physically distributed models to be location specific. There were still some aspects of a GIS that can enhance the linkage of the two models, which were not fully utilized in this study. As a powerful visualization tool, a GIS can be used to present the results of the MODFLOW model. Although results from the MODFLOW simulations were not transferred into a GIS, future research can help determine the best way to view MODFLOW?s results, by either grid (raster) or vector (polygon) coverages. In addition, new GIS software, such as Arc/Info 8.0, is Visual Basic compatible and would therefore be able to interact directly with the surface water and groundwater models. Consequently, the potential exists to automate some of the 111 processes that were accomplished manually for this project. These processes include determining river reach lengths, as well as the initial river discharge in the surface water model. 7.3.2 FORTRAN COMPILER This project also pointed to the importance of consistently using the same FORTRAN compiler if changes to the MODFLOW code are needed. Because of the possible incompatibility of the MODFLOW input files that would result by using a different FORTRAN compiler, the process of running MODFLOW was not fully automated. 7.4 Preliminary Results The results of the preliminary investigation into the effect of time-step in modeling groundwater and surface water interaction are presented in this section. As expected the smaller time-steps were able to capture the large changes in surface water discharge with more accuracy. Therefore, the fluctuations in the river and aquifer systems were larger when using a small time-step. The aquifer?s responses to the flood wave in Layers 1 and 2 were different. In Layer 1, aquifer heads responded more slowly to the flood wave with greater distance from the river. Furthermore, the fluctuations in the aquifer head farther from the river diminished in size. Aquifer heads in Layer 2, however, did not display either a reduction in fluctuations or a lag time in responding to the flood wave as a function of distance from the river. Also, aquifer head fluctuations in Layer 2 were much smaller than those in Layer 1. 112 A possible trend was discerned between the amount of water exchanged and time-step when examining the first 8, 16 and 24 days of the simulation. The smallest time-steps showed the smallest amount of river leakage into the aquifer. However, the dynamically linked models need to be validated. Furthermore, with the limited number of short simulations that were conducted, it is difficult to determine whether this trend can be applicable over a longer period of time and for other flood scenarios. Finally, an analysis of the lateral inflow and river characteristics indicated that the diffusive wave equation should be used in simulations that do not include the influence of a flood wave. Therefore, any future research may wish to change the surface water model to accommodate more diverse simulations. 7.5 Recommendations Future research of the Colorado River and the Carrizo-Wilcox aquifer may wish to consider the following recommendations: ? Modifying MODFLOW's FORTRAN code and compiling using a Lahey compiler to fully automate the dynamic linkage ? Conducting a sensitivity analysis of the results to the riverbed conductance used in MODFLOW's River Package and to the riverbed slopes used in Manning's equation ? Conducting simulations that exceed 32 days to assess whether the observed trend in river leakage rate is generally applicable 113 ? Considering a losing stream initial condition in which the river leakage caused by a flood wave would not return to the river as baseflow in the time-step following the peak discharge ? Introducing a well field near the river to assess the possible decrease in streamflow caused by well pumpage ? Using the diffusive wave approximation because of its applicability under non-flood scenarios 114 Appendix A: Surface Water Model Data 115 Table A1.River Discharge in Nodes Resulting from GIS Analysis and Model Calibration GIS Surface Water Dynamically Linked Node Analysis Model Calibration Model Calibration 0 65.70 77.30 77.30 1 65.75 77.35 77.35 2 65.82 77.42 77.42 3 65.86 77.46 77.46 4 66.38 77.98 77.98 5 66.39 77.99 77.99 6 66.41 78.01 78.01 7 66.46 78.06 78.06 8 66.58 78.18 78.18 9 66.59 78.19 78.19 10 66.61 78.21 78.21 11 66.62 78.22 78.22 12 66.63 78.23 78.23 13 66.64 78.24 78.24 14 66.79 78.39 78.40 15 66.80 78.40 78.40 16 67.16 78.76 78.76 17 67.31 78.91 78.91 18 67.35 78.95 78.95 19 68.01 79.61 79.61 20 68.04 79.64 79.65 21 68.38 79.98 79.99 22 68.41 80.01 80.01 23 68.47 80.07 80.07 24 68.49 80.09 80.09 25 68.53 80.13 80.14 26 68.54 80.14 80.14 27 68.61 80.21 80.22 28 68.64 80.24 80.25 29 69.77 81.37 81.37 30 69.79 81.39 81.39 31 69.89 81.49 81.49 32 70.02 81.62 81.62 33 70.03 81.63 81.64 34 70.05 81.65 81.66 35 71.88 83.48 83.49 River Discharge (m3/s) 116 Table A2. River Parameters in Par Worksheet Node Reach Length Distance from Inlet Reach Width Slope a (m) (m) (m) 0 0 0 76.2 1.71E-04 11.06 1 1120 1120 76.2 1.71E-04 11.06 2 2665 3785 76.2 1.19E-04 12.34 3 2445 6230 76.2 6.31E-05 14.93 4 3280 9510 76.2 1.28E-04 12.08 5 1890 11400 76.2 7.39E-05 14.24 6 1720 13120 76.2 7.76E-05 14.03 7 2260 15380 76.2 1.26E-04 12.12 8 2010 17390 76.2 6.95E-05 14.50 9 1990 19380 76.2 1.69E-04 11.10 10 1820 21200 76.2 3.84E-05 17.33 11 890 22090 76.2 7.84E-05 13.98 12 1120 23210 76.2 1.25E-04 12.17 13 1920 25130 76.2 1.03E-04 12.89 14 2550 27680 76.2 5.47E-05 15.58 15 710 28390 76.2 1.97E-04 10.61 16 1440 29830 76.2 9.70E-05 13.12 17 1860 31690 76.2 1.50E-04 11.51 18 1810 33500 76.2 1.54E-04 11.42 19 2540 36040 76.2 5.50E-05 15.56 20 1860 37900 76.2 7.51E-05 14.17 21 2490 40390 76.2 5.61E-05 15.47 22 1040 41430 76.2 1.34E-04 11.90 23 2450 43880 76.2 3.21E-05 18.29 24 1730 45610 76.2 4.54E-05 16.48 25 2600 48210 76.2 1.41E-04 11.72 26 970 49180 76.2 1.80E-04 10.90 27 2150 51330 76.2 6.50E-05 14.80 28 1930 53260 76.2 1.72E-04 11.05 29 3100 56360 76.2 8.12E-05 13.84 30 1870 58230 76.2 2.24E-04 10.21 31 6090 64320 76.2 6.31E-05 14.93 32 2210 66530 76.2 6.32E-05 14.92 33 1740 68270 76.2 8.03E-05 13.89 34 1640 69910 76.2 1.70E-04 11.08 35 4470 74380 76.2 6.25E-05 14.97 117 Table A3. Initial Conditions for Simulations Nodes River Discharge River Heads (m 3/s) (ft) (ft3/d) (m 3/s) 0 77.30 86.26 0 0 1 77.35 86.07 -628.76 -2.06E-04 2 77.42 86.19 -497.97 -1.63E-04 3 77.46 86.67 -450.86 -1.48E-04 4 77.98 84.53 -449.87 -1.47E-04 5 77.99 84.42 -366.33 -1.20E-04 6 78.01 83.84 -377.56 -1.24E-04 7 78.06 82.28 -483.69 -1.59E-04 8 78.18 82.76 -402.72 -1.32E-04 9 78.19 80.29 -544.30 -1.78E-04 10 78.21 82.86 -439.58 -1.44E-04 11 78.22 80.66 -457.05 -1.50E-04 12 78.23 79.36 -380.79 -1.25E-04 13 78.24 79.33 -466.57 -1.53E-04 14 78.40 80.28 -375.20 -1.23E-04 15 78.40 76.89 -284.21 -9.31E-05 16 78.76 77.93 -536.52 -1.76E-04 17 78.91 76.53 -706.95 -2.32E-04 18 78.95 75.56 -323.84 -1.06E-04 19 79.61 77.14 -351.18 -1.15E-04 20 79.65 75.86 -445.48 -1.46E-04 21 79.99 76.20 -364.60 -1.19E-04 22 80.01 73.61 -716.33 -2.35E-04 23 80.07 76.97 88.16 2.89E-05 24 80.09 75.63 -511.04 -1.67E-04 25 80.14 72.54 -354.88 -1.16E-04 26 80.14 70.84 -758.11 -2.48E-04 27 80.22 72.60 -174.69 -5.73E-05 28 80.25 69.91 -821.38 -2.69E-04 29 81.37 70.56 -1179.61 -3.87E-04 30 81.39 67.54 -798.66 -2.62E-04 31 81.49 69.02 -750.85 -2.46E-04 32 81.62 67.76 -1414.47 -4.64E-04 33 81.64 66.68 -277.97 -9.11E-05 34 81.66 64.53 -1349.15 -4.42E-04 35 83.49 66.08 -2863.41 -9.38E-04 River Leakage Rates 118 Table A4. Nodes and Corresponding River Package Cells in MODFLOW Nodes Layer Row Column 0 - - - 1 1 21 5 2 1 21 6 3 1 22 7 4 1 22 8 5 1 23 8 6 1 24 8 7 1 25 8 8 1 26 9 9 1 26 10 10 1 27 11 11 1 28 11 12 1 28 10 13 1 29 10 14 1 30 9 15 1 31 9 16 1 31 8 17 1 32 8 18 1 32 9 19 1 32 10 20 1 33 10 21 1 34 10 22 1 35 10 23 1 35 11 24 1 35 12 25 1 36 12 26 1 36 11 27 1 37 12 28 1 37 13 29 1 38 14 30 1 39 13 31 1 39 11,12* 32 1 40 13 33 1 40 12 34 1 40 11 35 1 41 11 River Package Cells * Two River Package cells were modeled as one reach in the surface water model 119 Table A5. Inlet River Discharge for Various Time-Steps Simulation Day Date 1-Day* 2-Day 4-Day 8-Day 16-Day 0 1/8/1991 11.19 14.17 412.01 266.39 227.05 1 1/9/1991 17.16 14.17 412.01 266.39 227.05 2 1/10/1991 586.15 809.85 412.01 266.39 227.05 3 1/11/1991 1033.56 809.85 412.01 266.39 227.05 4 1/12/1991 223.13 146.82 120.77 266.39 227.05 5 1/13/1991 70.51 146.82 120.77 266.39 227.05 6 1/14/1991 51.25 94.72 120.77 266.39 227.05 7 1/15/1991 138.18 94.72 120.77 266.39 227.05 8 1/16/1991 171.88 125.73 275.66 187.70 227.05 9 1/17/1991 79.57 125.73 275.66 187.70 227.05 10 1/18/1991 267.87 425.60 275.66 187.70 227.05 11 1/19/1991 583.32 425.60 275.66 187.70 227.05 12 1/20/1991 230.78 153.05 99.75 187.70 227.05 13 1/21/1991 75.32 153.05 99.75 187.70 227.05 14 1/22/1991 52.10 46.44 99.75 187.70 227.05 15 1/23/1991 40.78 46.44 99.75 187.70 227.05 16 1/24/1991 43.61 49.13 46.72 36.66 73.76 17 1/25/1991 54.65 49.13 46.72 36.66 73.76 18 1/26/1991 52.39 44.32 46.72 36.66 73.76 19 1/27/1991 36.25 44.32 46.72 36.66 73.76 20 1/28/1991 30.02 29.17 26.60 36.66 73.76 21 1/29/1991 28.32 29.17 26.60 36.66 73.76 22 1/30/1991 24.98 24.04 26.60 36.66 73.76 23 1/31/1991 23.11 24.04 26.60 36.66 73.76 24 2/1/1991 21.95 21.68 54.42 110.85 73.76 25 2/2/1991 21.41 21.68 54.42 110.85 73.76 26 2/3/1991 19.99 87.16 54.42 110.85 73.76 27 2/4/1991 154.33 87.16 54.42 110.85 73.76 28 2/5/1991 427.58 286.85 167.28 110.85 73.76 29 2/6/1991 146.11 286.85 167.28 110.85 73.76 30 2/7/1991 56.63 47.71 167.28 110.85 73.76 31 2/8/1991 38.79 47.71 167.28 110.85 73.76 * 1-day average is equivalent to the original data provided by USGS gage 8160400 120 Appendix B: Interface Visual Basic Code 121 Procedure cmdExit_Click Private Sub cmdExit_Click() 'Function: Exits the model 'Procedure runs when user click Exit button on frmInput form End End Sub Procedure cmdOK_Click Private Sub cmdOK_Click() 'Function: Performs dynamic linkage of MODFLOW and surface water models 'Procedure runs when user clicks OK on frmInput form Dim iTimestep As Integer Dim sPath As String Dim iLoopcount As Integer Dim iTotTime As Double Dim dMultiFactor As Double Dim iTime As Integer Application.DisplayAlerts = False iLoopcount = 0 iTotTime = txtTotTime.Text iTimestep = txtTimestep.Text sPath = txtPathname.Text iTime = 0 Initialize sPath, iLoopcount Do iLoopcount = iLoopcount + 1 iTime = iTime + iTimestep Outputsheet iTimestep, iTime, iLoopcount Rename sPath, iLoopcount Open_Bas iTimestep, sPath, iLoopcount Open_Riv sPath Run_MODFLOW sPath Output_qa sPath, iLoopcount Loop Until iTime >= iTotTime Application.DisplayAlerts = True End End Sub 122 Sub-Procedure Initialize 'Function: Initializes surface water model and the three output worksheets for new 'simulation Sub Initialize(sPath As String, iLoopcount As Integer) 'Copies initial leakage rates from Par worksheet to qa,h_fix worksheet Sheets("Par").Select Range("A28:AJ28").Select Selection.Copy Sheets("qa,h_fix").Select Range("D5").Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=True 'Copies initial flowrate from Par worksheet to Q worksheet Sheets("Par").Select Range("A16:AJ16").Select Selection.Copy Sheets("Q").Select Range("D3").Select ActiveSheet.Paste Windows("kinematic8").Activate 'Clears flowrate worksheet of any results from previous simulation Sheets("flowrate").Select Rows("4:100").Select Selection.ClearContents 'Clears riv_head worksheet of any results from previous simulation Sheets("riv_head").Select Rows("4:100").Select Selection.ClearContents 'Clears aq_head worksheet of any results from previous simulation Sheets("aq_head").Select Rows("4:100").Select Selection.ClearContents 123 Sub-Procedure Surface_Run 'Function: Sets new inlet river discharges as the upper boundary condition for a 'one-day simulation Private Sub Surface_Run(n As Integer, iTime As Integer, iTimestep As Integer) Dim iDay As Integer iDay = n + iTime - iTimestep If iTimestep = 2 Then Sheets("inlet2").Select ElseIf iTimestep = 4 Then Sheets("inlet4").Select ElseIf iTimestep = 8 Then Sheets("inlet8").Select ElseIf iTimestep = 16 Then Sheets("inlet16").Select End If Range(Cells(3, iDay + 4), Cells(483, iDay + 4)).Select Selection.Copy Sheets("Q").Select Range("D3").Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False End Sub 124 Sub-Procedure Outputsheet 'Function: Runs Surface_Run for the number of days in the interface time-step (iTimestep). Also. copies the river heads and discharge resulting from the simulation into the flowrate and riv_head worksheets, respectively. Private Sub Outputsheet (iTimestep As Integer, iTime As Integer, iLoopcount As Integer) Dim n As Integer Dim yo As Double Dim iCol As Integer 'Initialize parameters n = 0 Do If iTimestep = 2 Or iTimestep = 4 Or iTimestep = 8 Or iTimestep = 16 Then Surface_Run n, iTime, iTimestep End If Sheets("Q").Select Range("D483:AN483").Select Selection.Copy ActiveWindow.LargeScroll Down:=-10 Range("D3").Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False n = n + 1 Loop Until n >= iTimestep 'Copies river flowrate resulting from simulation from the Q worksheet to the flowrate sheet Range("D482:AM482").Copy Sheets("flowrate").Select Cells(iLoopcount + 3, 3).Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False 'Copies river heads resulting from simulation from the h worksheet into qa,h_fix worksheet Sheets("h").Select Range("E482:AM482").Copy yo = Cells(482, 4).Value Sheets("qa,h_fix").Select Range("M5").Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=True 'Copies river heads for nodes 1-36 from sheet qa,h_fix, and pastes in riv_head sheet Sheets("qa,h_fix").Select Range("P5:P39").Select Selection.Copy Sheets("riv_head").Select Range(Cells(iLoopcount + 3, 4), Cells(iLoopcount + 3, 38)).Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=True 125 Sub-Procedure Outputsheet cont'd 'Writes the river hydraulic head for reach 0 in the river head reach 0 column of riv_head sheet Cells(iLoopcount + 3, 3).Value = yo * 3.281 + 79.8 'Writes the time total simulation time in the second column of flowrate the sheet Sheets("flowrate").Select Cells(iLoopcount + 3, 2).Value = iTime 'Writes the loop index number in the first column of flowrate sheet Cells(iLoopcount + 3, 1).Value = iLoopcount 'Writes the total simulation time in the second column of the riv_head sheet Sheets("riv_head").Select Cells(iLoopcount + 3, 2).Value = iTime 'Writes the loop index number in the first column of riv_head sheet Cells(iLoopcount + 3, 1).Value = iLoopcount 'Writes the total simulation time 'in the second column of the aq_head sheet Sheets("aq_head").Select Cells(iLoopcount + 3, 2).Value = iTime 'Writes the loop index number in the first column of aq_head sheet Cells(iLoopcount + 3, 1).Value = iLoopcount End Sub 126 Sub-Procedure Rename 'Function: Initializes the starting heads for the first MODFLOW run (i.e. iLoopcount = 1). Also, 'renames *.hed file at the end of MODFLOW run to the *.shd (starting head file). The *.shd file 'contains the starting heads used in the beginning of the next MODFLOW run. Private Sub Rename(sPath As String, iLoopcount As Integer) If iLoopcount = 1 Then FileCopy sPath & "\modflow\initial.bas", sPath & "\modflow\Cw409_s2.bas" FileCopy sPath & "\modflow\initial.nam", sPath & "\modflow\Cw409_s2.nam" ElseIf iLoopcount = 2 Then FileCopy sPath & "\modflow\standard.bas", sPath & "\modflow\Cw409_s2.bas" FileCopy sPath & "\modflow\standard.nam", sPath & "\modflow\Cw409_s2.nam" FileCopy sPath & "\modflow\Cw409_s2.hed", sPath & "\modflow\Cw409_s2.shd" ElseIf iLoopcount > 2 Then FileCopy sPath & "\modflow\Cw409_s2.hed", sPath & "\modflow\Cw409_s2.shd" End If End Sub 127 Sub-Procedure Open_Bas Sub Open_Bas(iTimestep As Integer, sPath As String, iLoopcount As Integer) ' Function: Enters the time-step for the MODFLOW run into the Basic Package 'It's important that the Basic Package initially be saved as a *.csv (comma delimited) file Workbooks.OpenText Filename:=sPath & "\modflow\CW409_S2.BAS", Origin _ :=xlWindows, StartRow:=1, DataType:=xlDelimited, TextQualifier:= _ xlDoubleQuote, ConsecutiveDelimiter:=True, Tab:=False, Semicolon:=False, _ Comma:=True, Space:=False, Other:=False, FieldInfo:=Array(Array(1, 1), _ Array(2, 1), Array(3, 1), Array(4, 1), Array(5, 1), Array(6, 1), Array(7, 1), Array(8, 1), _ Array(9, 1)) 'Sets MODFLOW variable NPER (number of stress periods) = 1 Range("E3").Select ActiveCell.FormulaR1C1 = "1" 'Sets MODFLOW variable PERLEN (length of stress period) equal to user specified time-step If iLoopcount = 1 Then Range("A437").Select ActiveCell.Value = iTimestep Else Range("A227").Select ActiveCell.Value = iTimestep End If 'Saves and closes Basic Package ActiveWorkbook.SaveAs Filename:=sPath & "\modflow\CW409_S2.BAS", _ FileFormat:=xlCSV, CreateBackup:=False ActiveWorkbook.Close End Sub 128 Sub-Procedure Open_Riv 'Funtion: Opens River Package and inserts new river heads from sheet qa,h_fix 'River Package (*.riv) must be saved as a csv (comma delimited) file Sub Open_Riv(sPath As String) Workbooks.OpenText Filename:=sPath & "\modflow\Cw409_S2.riv", Origin:= _ xlWindows, StartRow:=1, DataType:=xlDelimited, TextQualifier:= _ xlDoubleQuote, ConsecutiveDelimiter:=True, Tab:=False, Semicolon:=False, _ Comma:=True, Space:=False, Other:=False, FieldInfo:=Array(Array(1, 1), _ Array(2, 1), Array(3, 1), Array(4, 1), Array(5, 1), Array(6, 1), Array(7, 1)) Windows("kinematic8.xls").Activate Sheets("qa,h_fix" ).Select 'Copies cells which contain the new river heads from qa,h_fix worksheet Range("T5:T40").Select Selection.Copy Windows("CW409_S2.RIV").Activate 'Pastes onto cell D67 of the River Package Range("D67").Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False Range("D67:D108").Select Selection.NumberFormat = "0.00" 'Saves as a .csv file and closes ActiveWorkbook.SaveAs Filename:=sPath & "\modflow\CW409_S2.RIV", FileFormat _ :=xlCSV, CreateBackup:=False ActiveWorkbook.Close End Sub Sub-Procedure Run_MODFLOW 'Function: Runs MODFLOW. This sub-procedure is not fully automated and requires user to 'enter in the MODFLOW Name file (cw409_S2.nam). Private Sub Run_MODFLOW(sPath As String) Dim Start As Double Dim ID As Integer ID = Shell(sPath & "\modflow\Modflw96.exe") MsgBox ("RUN MODFLOW") End Sub 129 Sub-Procedure Output_qa 'Function: Transfers river leakage rates from MODFLOW listing file (*.out) to qa,h_fix sheet 'Also transfers aquifer heads in Cross Sections 1 and 2 to aq_head sheet Private Sub Output_qa(sPath As String, iLoopcount As Integer) Dim iRow As Integer Dim iCount As Integer Dim iRowaqhead As Integer Dim lay1 As Boolean 'initialize iRow iRow = 0 'Opens listing file (.out) Workbooks.OpenText Filename:= _ sPath & "\modflow\Cw409_s2.out", Origin:=xlWindows, _ StartRow:=1, DataType:=xlDelimited, TextQualifier:=xlDoubleQuote, _ ConsecutiveDelimiter:=True, Tab:=False, Semicolon:=False, Comma:=False, _ Space:=True, Other:=False, FieldInfo:=Array(Array(1, 1), Array(2, 1), Array(3 _ , 1), Array(4, 1), Array(5, 1), Array(6, 1), Array(7, 1), Array(8, 1)) 'Uses Excel's Find command to locate Reach 1 (or REACH 65 in revised BEG model) and copies river leakage rates. Searches for the word " leakage" in column C and assigns that row as iRow Columns("C:C").Select iRow = Selection.Find(What:="leakage", After:=ActiveCell, LookIn:=xlFormulas, _ LookAt:=xlPart, SearchOrder:=xlByRows, SearchDirection:=xlNext, _ MatchCase:=False).Row Range(Cells(iRow + 65, 11), Cells(iRow + 100, 11)).Select Selection.Copy Windows("kinematic8.xls").Activate ' Pastes river leakage rates into sheet qa,h_fix Sheets("qa,h_fix").Select Range("D5").Select ActiveSheet.Paste ' Copies river leaksge rates for reaches 14, 27, 28 and for the all of the 35 reaches being modeled and pastes into aq_head sheet 'Reach 14 Sheets("qa,h_fix").Select Range("D16").Select Selection.Copy Sheets("aq_head").Select ' Range("AO3").Select Cells(iLoopcount + 3, 41).Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False 'Reach 27 Sheets("qa,h_fix").Select Range("D34").Select 130 Selection.Copy Sheets("aq_head").Select 'Range("AQ3").Select Cells(iLoopcount + 3, 42).Select ActiveSheet.Paste 'Reach 28 Sheets("qa,h_fix").Select Range("D39").Select Selection.Copy Sheets("aq_head").Select ' Range("AR3").Select Cells(iLoopcount + 3, 43).Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False 'All 35 reaches Sheets("qa,h_fix").Select Range("D41").Select Selection.Copy Sheets("aq_head").Select ' Range("AR3").Select Cells(iLoopcount + 3, 44).Select Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _ False, Transpose:=False ' Locates the final aquifer heads in Layer 1 resulting from the simulation in the Listing file. 'Designates iRow as the row where "end" appears in Column G. Windows("Cw409_s2.out").Activate iRow = 0 Columns("G:G").Select iRow = Selection.Find(What:="END", After:=ActiveCell, LookIn:=xlFormulas, _ LookAt:=xlPart, SearchOrder:=xlByRows, SearchDirection:=xlNext, _ MatchCase:=False).Row iRowaqhead = iRow + 193 ' Locates aquifer heads in Cross Sections 1 and 2 by referencing iRow lay1 = True makegrid iRowaqhead, iLoopcount outputhead iRowaqhead, iLoopcount, lay1 iRowaqhead = iRowaqhead + 216 lay1 = False makegrid iRowaqhead, iLoopcount outputhead iRowaqhead, iLoopcount, lay1 ' Closes listing file Windows("CW409_S2.OUT").Activate ActiveWorkbook.Close End Sub 131 Sub-Procedure Makegrid Function: Reformats aquifer heads in the MODFLOW listing file to resemble the domain grid. Private Sub makegrid(iRowaqhead As Integer, iLoopcount As Integer) Dim iRow2 As Integer Dim Count As Integer Dim iFinish As Integer Dim sFile As String If iLoopcount = 0 Then sFile = "initial.out" Else sFile = "Cw409_s2.out" End If Windows(sFile).Activate iRow2 = iRowaqhead iFinish = iRow2 + 22 Do Range(Cells(iRow2 + 1, 2), Cells(iRow2 + 1, 10)).Select Selection.Cut Cells(iRow2, 12).Select ActiveSheet.Paste Range(Cells(iRow2 + 2, 2), Cells(iRow2 + 2, 10)).Select Selection.Cut Cells(iRow2, 21).Select ActiveSheet.Paste Range(Cells(iRow2 + 3, 2), Cells(iRow2 + 3, 10)).Select Selection.Cut Cells(iRow2, 30).Select ActiveSheet.Paste Range(Cells(iRow2 + 4, 2), Cells(iRow2 + 4, 10)).Select Selection.Cut Cells(iRow2, 39).Select ActiveSheet.Paste Range(Cells(iRow2 + 5, 2), Cells(iRow2 + 5, 10)).Select Selection.Cut Cells(iRow2, 48).Select ActiveSheet.Paste Range(Cells(iRow2 + 6, 2), Cells(iRow2 + 6, 10)).Select Selection.Cut Cells(iRow2, 57).Select ActiveSheet.Paste Range(Cells(iRow2 + 7, 2), Cells(iRow2 + 7, 10)).Select Selection.Cut Cells(iRow2, 66).Select ActiveSheet.Paste Range(Cells(iRow2 + 8, 2), Cells(iRow2 + 8, 10)).Select 132 Sub-Procedure Makegrid cont'd Selection.Cut Cells(iRow2, 75).Select ActiveSheet.Paste Count = 0 Do Count = Count + 1 Rows(iRow2 + 1).Delete Loop Until (Count = 8) iRow2 = iRow2 + 1 Loop Until iRow2 = iFinish End Sub 133 Sub-Procedure Outputhead 'Function: Extracts aquifer heads from Cross Section 1 and 2 from the MODFLOW listing file Sub outputhead(iRow As Integer, iLoopcount As Integer, lay1 As Boolean) Dim sFile As String Dim iRow2 As Integer iRow2 = iRow sFile = "Cw409_s2.out" 'Cross Section 1 Windows(sFile).Activate Range(Cells(iRow2 + 9, 7), Cells(iRow2 + 9, 15)).Select Selection.Copy Windows("kinematic8.xls").Activate Sheets("aq_head").Select If lay1 = True Then Cells(iLoopcount + 3, 3).Select Else Cells(iLoopcount + 3, 12).Select End If ActiveSheet.Paste 'Cross Section 2 Windows(sFile).Activate Range(Cells(iRow2 + 16, 10), Cells(iRow2 + 16, 19)).Select Selection.Copy Windows("kinematic8.xls").Activate If lay1 = True Then Cells(iLoopcount + 3, 21).Select Else Cells(iLoopcount + 3, 31).Select End If ActiveSheet.Paste End Sub 134 Appendix C: Simulation Results 135 Table C1. Cross Section 1 Aquifer Heads in Layer 1 from 2-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 1 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 110.628 106.773 99.929 88.751 80.288 85.421 91.612 97.677 103.902 2 110.627 106.762 99.813 87.728 74.555 84.360 91.415 97.643 103.896 4 110.633 106.806 100.340 93.106 108.324 89.893 92.310 97.777 103.917 6 110.643 106.865 100.623 92.312 84.573 89.256 92.783 97.945 103.957 8 110.654 106.919 100.720 91.307 81.475 88.293 92.932 98.083 104.005 10 110.667 106.964 100.767 91.028 83.353 88.011 92.995 98.190 104.056 12 110.683 107.029 101.055 93.200 96.448 90.254 93.466 98.348 104.118 14 110.702 107.091 101.193 92.540 84.889 89.648 93.677 98.495 104.185 16 110.719 107.132 101.116 90.902 77.926 87.971 93.523 98.572 104.244 18 110.733 107.150 100.957 89.900 78.141 86.879 93.234 98.582 104.288 20 110.744 107.148 100.767 89.171 77.722 86.052 92.896 98.538 104.314 22 110.751 107.129 100.549 88.427 76.295 85.210 92.518 98.450 104.322 24 110.754 107.094 100.324 87.827 75.743 84.518 92.136 98.330 104.311 26 110.753 107.048 100.106 87.366 75.473 83.976 91.774 98.189 104.284 28 110.750 107.007 100.019 88.021 80.946 84.600 91.637 98.073 104.250 30 110.748 106.994 100.199 90.235 91.117 86.889 91.957 98.049 104.225 32 110.746 106.984 100.217 89.313 78.027 86.016 91.997 98.034 104.206 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction. 136 Table C2. Cross Section 1 Aquifer Heads in Layer 2 from 2-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 2 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 101.547 101.651 101.765 101.890 102.039 102.170 102.238 102.256 102.355 2 101.527 101.630 101.742 101.866 102.012 102.144 102.213 102.232 102.334 4 101.635 101.745 101.864 101.999 102.156 102.284 102.349 102.360 102.450 6 101.635 101.743 101.859 101.988 102.139 102.270 102.336 102.352 102.445 8 101.627 101.734 101.849 101.977 102.126 102.258 102.325 102.341 102.436 10 101.628 101.735 101.851 101.979 102.129 102.260 102.327 102.343 102.438 12 101.678 101.787 101.906 102.039 102.193 102.323 102.389 102.401 102.491 14 101.676 101.784 101.901 102.031 102.182 102.313 102.380 102.394 102.487 16 101.651 101.757 101.873 102.000 102.149 102.281 102.349 102.365 102.460 18 101.632 101.738 101.853 101.980 102.128 102.260 102.328 102.345 102.442 20 101.614 101.720 101.834 101.960 102.109 102.241 102.309 102.327 102.424 22 101.594 101.699 101.813 101.938 102.086 102.218 102.286 102.305 102.404 24 101.575 101.679 101.792 101.917 102.064 102.196 102.265 102.284 102.384 26 101.557 101.661 101.774 101.898 102.045 102.177 102.245 102.265 102.366 28 101.561 101.666 101.779 101.905 102.054 102.185 102.253 102.272 102.372 30 101.600 101.707 101.823 101.952 102.104 102.234 102.301 102.317 102.413 32 101.586 101.691 101.805 101.931 102.079 102.211 102.279 102.297 102.396 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction. 137 Table C3. Cross Section 2 Aquifer Heads in Layer 1 from 2-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 95.616 90.014 83.018 76.919 72.604 69.927 79.385 87.997 94.995 100.670 2 95.612 89.992 82.882 76.094 67.298 65.952 78.659 87.888 94.980 100.668 4 95.624 90.077 83.505 80.444 99.063 89.691 82.507 88.393 95.040 100.675 6 95.651 90.197 83.880 80.062 76.715 72.993 82.124 88.691 95.124 100.692 8 95.688 90.310 84.042 79.383 73.740 70.771 81.494 88.813 95.202 100.714 10 95.730 90.408 84.140 79.207 75.503 72.081 81.324 88.883 95.269 100.739 12 95.783 90.543 84.520 81.021 87.861 81.319 82.915 89.183 95.362 100.770 14 95.844 90.678 84.735 80.637 76.987 73.200 82.547 89.347 95.454 100.806 16 95.903 90.772 84.694 79.379 70.424 68.289 81.419 89.303 95.518 100.841 18 95.956 90.820 84.538 78.532 70.598 68.422 80.671 89.170 95.549 100.872 20 95.997 90.827 84.327 77.866 70.211 68.133 80.093 88.996 95.552 100.896 22 96.024 90.796 84.069 77.172 68.877 67.137 79.496 88.788 95.529 100.912 24 96.037 90.734 83.789 76.579 68.373 66.763 78.994 88.565 95.485 100.920 26 96.034 90.647 83.508 76.094 68.121 66.572 78.587 88.344 95.425 100.920 28 96.022 90.566 83.373 76.513 73.216 70.376 78.978 88.242 95.369 100.913 30 96.011 90.538 83.561 78.276 82.810 77.551 80.547 88.401 95.352 100.908 32 96.001 90.519 83.575 77.609 70.531 68.378 79.949 88.417 95.341 100.904 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 138 Table C4. Cross Section 2 Aquifer Heads in Layer 2 from 2-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 98.268 98.074 97.901 97.677 97.583 97.596 97.676 97.789 97.873 97.978 2 98.246 98.050 97.874 97.648 97.552 97.572 97.659 97.779 97.865 97.971 4 98.365 98.182 98.024 97.816 97.730 97.714 97.761 97.844 97.915 98.018 6 98.350 98.159 97.990 97.767 97.665 97.665 97.731 97.830 97.906 98.010 8 98.335 98.141 97.969 97.743 97.641 97.645 97.717 97.821 97.900 98.004 10 98.336 98.142 97.971 97.745 97.644 97.648 97.718 97.822 97.900 98.005 12 98.389 98.201 98.036 97.817 97.719 97.708 97.762 97.851 97.923 98.026 14 98.379 98.187 98.017 97.790 97.685 97.682 97.746 97.842 97.917 98.021 16 98.349 98.154 97.979 97.748 97.642 97.647 97.721 97.827 97.906 98.010 18 98.332 98.136 97.961 97.731 97.627 97.635 97.711 97.819 97.900 98.004 20 98.318 98.121 97.946 97.716 97.614 97.624 97.702 97.813 97.895 97.999 22 98.300 98.104 97.928 97.698 97.597 97.610 97.691 97.805 97.888 97.993 24 98.285 98.088 97.911 97.682 97.583 97.599 97.682 97.798 97.883 97.988 26 98.271 98.073 97.897 97.668 97.571 97.589 97.674 97.792 97.878 97.983 28 98.280 98.085 97.912 97.686 97.592 97.605 97.684 97.798 97.882 97.987 30 98.323 98.133 97.965 97.746 97.653 97.654 97.720 97.821 97.899 98.003 32 98.302 98.106 97.932 97.705 97.605 97.617 97.696 97.807 97.890 97.995 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 139 Table C5. Leakage Rates from 2-Day Time-Step Simulation Leakage Rate (ft3/d) Time (days) Reach14 Reach 27 Reach 28 Total 0 -375.20 -174.69 -821.38 -21216.42 2 -1182.82 -1842.06 -2028.62 -237321.32 4 4145.75 9477.00 6198.37 1224020.27 6 -2549.99 -5508.15 -4947.08 -757116.01 8 -722.14 -966.11 -1506.98 -206757.09 10 -137.39 339.90 -527.27 28122.22 12 1537.29 3797.98 1964.95 478992.32 14 -1393.40 -2724.72 -2871.44 -388353.93 16 -1263.92 -2114.46 -2349.31 -317512.74 18 -527.87 -402.69 -1030.72 -49928.88 20 -576.54 -529.55 -1099.28 -45093.05 22 -714.59 -822.76 -1289.16 -82071.65 24 -633.40 -621.18 -1114.81 -53807.83 26 -599.79 -546.89 -1044.50 -40479.20 28 212.32 1142.19 187.53 177597.38 30 1109.75 2897.35 1408.93 383904.14 32 -1760.59 -3384.36 -3237.97 -460118.96 140 Table C6a. River Discharge at Finite Difference Nodes for a 2-Day Time-Step Simulation River Discharge (m3/s) Time (days) 0 1 2 3 4 5 6 7 8 0 77.300 77.348 77.416 77.456 77.982 77.991 78.008 78.062 78.181 2 14.172 14.221 14.289 14.329 14.856 14.864 14.881 14.935 15.054 4 809.854 809.908 809.982 810.033 810.568 810.589 810.618 810.673 810.794 6 146.821 146.836 146.875 146.848 147.329 147.268 147.213 147.261 147.372 8 94.719 94.788 94.875 94.954 95.508 95.556 95.614 95.672 95.797 10 125.726 125.779 125.851 125.902 126.436 126.457 126.487 126.542 126.661 12 425.598 425.645 425.712 425.749 426.274 426.281 426.296 426.349 426.468 14 153.051 153.086 153.143 153.155 153.664 153.645 153.633 153.684 153.800 16 46.439 46.498 46.575 46.635 47.175 47.204 47.242 47.298 47.420 18 49.129 49.185 49.260 49.317 49.855 49.881 49.917 49.972 50.092 20 44.315 44.364 44.433 44.474 45.002 45.013 45.032 45.086 45.205 22 29.166 29.215 29.283 29.324 29.852 29.862 29.880 29.934 30.053 24 24.041 24.091 24.160 24.203 24.732 24.744 24.764 24.818 24.938 26 21.676 21.725 21.794 21.836 22.363 22.374 22.393 22.447 22.566 28 87.158 87.207 87.276 87.316 87.844 87.853 87.871 87.925 88.044 30 286.847 286.890 286.953 286.982 287.501 287.499 287.504 287.557 287.675 32 47.713 47.751 47.810 47.827 48.339 48.325 48.318 48.370 48.487 141 Table C6b. River Discharge at Finite Difference Nodes for a 2-Day Time-Step Simulation River Discharge (m3/s)) Time (days) 9 10 11 12 13 14 15 16 17 0 78.192 78.206 78.220 78.226 78.240 78.393 78.396 78.761 78.908 2 15.066 15.080 15.094 15.101 15.116 15.269 15.272 15.637 15.785 4 810.806 810.821 810.835 810.844 810.858 811.012 811.015 811.380 811.528 6 147.377 147.390 147.403 147.402 147.415 147.567 147.569 147.933 148.079 8 95.813 95.827 95.842 95.854 95.869 96.023 96.026 96.392 96.540 10 126.674 126.688 126.702 126.710 126.724 126.878 126.881 127.246 127.394 12 426.479 426.493 426.507 426.514 426.528 426.681 426.684 427.049 427.197 14 153.809 153.823 153.837 153.840 153.854 154.007 154.009 154.374 154.521 16 47.433 47.448 47.462 47.471 47.486 47.639 47.643 48.008 48.156 18 50.105 50.120 50.134 50.143 50.158 50.311 50.314 50.680 50.827 20 45.217 45.231 45.245 45.253 45.267 45.420 45.423 45.788 45.936 22 30.065 30.079 30.094 30.101 30.115 30.268 30.271 30.637 30.784 24 24.950 24.964 24.978 24.986 25.000 25.153 25.156 25.522 25.670 26 22.578 22.592 22.606 22.614 22.628 22.781 22.784 23.149 23.297 28 88.056 88.070 88.085 88.092 88.106 88.259 88.262 88.628 88.775 30 287.685 287.699 287.713 287.719 287.733 287.886 287.889 288.254 288.401 32 48.496 48.510 48.524 48.528 48.542 48.695 48.698 49.063 49.210 142 Table C6c. River Discharge at Finite Difference Nodes for a 2-Day Time-Step Simulation River Discharge (m3/s) Time (days) 18 19 20 21 22 23 24 25 26 0 78.945 79.611 79.642 79.982 80.010 80.069 80.089 80.134 80.137 2 15.821 16.487 16.518 16.859 16.887 16.948 16.969 17.018 17.022 4 811.565 812.231 812.263 812.603 812.632 812.692 812.712 812.759 812.762 6 148.115 148.779 148.808 149.146 149.173 149.229 149.248 149.289 149.290 8 96.578 97.245 97.277 97.618 97.647 97.709 97.730 97.779 97.782 10 127.431 128.096 128.128 128.468 128.497 128.557 128.577 128.623 128.626 12 427.234 427.899 427.930 428.270 428.299 428.358 428.378 428.424 428.427 14 154.557 155.222 155.253 155.592 155.620 155.678 155.698 155.742 155.744 16 48.194 48.859 48.891 49.232 49.261 49.321 49.342 49.389 49.392 18 50.865 51.530 51.562 51.902 51.931 51.992 52.012 52.059 52.062 20 45.973 46.638 46.670 47.010 47.038 47.098 47.118 47.164 47.167 22 30.821 31.487 31.518 31.858 31.887 31.947 31.967 32.013 32.016 24 25.707 26.372 26.404 26.744 26.772 26.832 26.852 26.898 26.901 26 23.334 24.000 24.031 24.371 24.400 24.459 24.479 24.526 24.528 28 88.812 89.478 89.509 89.849 89.878 89.938 89.957 90.004 90.006 30 288.438 289.104 289.135 289.474 289.503 289.562 289.582 289.627 289.629 32 49.246 49.911 49.942 50.281 50.309 50.368 50.387 50.432 50.434 143 Table C6d. River Discharge at Finite Difference Nodes for a 2-Day Time-Step Simulation River Discharge (m3/s) Time (days) 27 28 29 30 31 32 33 34 35 0 80.211 80.242 81.366 81.386 81.485 81.615 81.632 81.647 83.479 2 17.101 17.138 18.284 18.317 18.644 18.883 18.997 19.095 21.414 4 812.836 812.868 813.994 814.014 814.117 814.248 814.266 814.281 816.113 6 149.362 149.391 150.511 150.528 150.611 150.741 150.754 150.766 152.598 8 97.858 97.891 99.018 99.040 99.150 99.281 99.301 99.317 101.149 10 128.700 128.732 129.857 129.877 129.979 130.109 130.127 130.142 131.974 12 428.500 428.532 429.656 429.676 429.775 429.906 429.923 429.937 431.769 14 155.817 155.847 156.970 156.989 157.082 157.212 157.228 157.241 159.073 16 49.467 49.499 50.625 50.646 50.751 50.882 50.901 50.916 52.751 18 52.137 52.169 53.294 53.315 53.418 53.549 53.567 53.582 55.414 20 47.241 47.273 48.397 48.417 48.518 48.648 48.666 48.680 50.513 22 32.090 32.121 33.246 33.266 33.370 33.502 33.521 33.537 35.384 24 26.976 27.007 28.132 28.152 28.258 28.390 28.410 28.427 30.277 26 24.603 24.634 25.759 25.779 25.884 26.016 26.036 26.052 27.899 28 90.080 90.112 91.237 91.257 91.358 91.488 91.505 91.520 93.352 30 289.703 289.734 290.858 290.877 290.975 291.106 291.122 291.137 292.969 32 50.507 50.538 51.661 51.680 51.775 51.906 51.922 51.936 53.770 144 Table C7a. River Hydraulic Heads at Finite Difference Nodes for a 2-Day Time-Step Simulation River Head (ft) Time (days) 0 1 2 3 4 5 6 7 8 0 86.262 86.070 86.194 86.670 84.526 84.420 83.838 82.284 82.756 2 82.137 81.942 81.591 81.105 80.050 79.145 78.643 77.797 77.395 4 106.278 106.079 108.520 113.670 106.361 110.153 109.190 104.192 108.961 6 89.304 89.105 89.578 90.759 87.827 88.306 87.664 85.589 86.708 8 87.106 86.909 87.131 87.806 85.445 85.505 84.910 83.210 83.863 10 88.459 88.262 88.639 89.627 86.915 87.236 86.613 84.681 85.623 12 97.798 97.600 99.058 102.227 97.103 99.241 98.440 94.902 97.847 14 89.544 89.345 89.847 91.086 88.093 88.622 87.977 85.860 87.032 16 84.564 84.368 84.296 84.376 82.677 82.242 81.695 80.432 80.543 18 84.728 84.531 84.478 84.596 82.854 82.451 81.900 80.609 80.754 20 84.432 84.235 84.147 84.195 82.530 82.068 81.522 80.283 80.364 22 83.404 83.207 83.001 82.810 81.416 80.755 80.228 79.165 79.028 24 83.009 82.813 82.562 82.279 80.989 80.253 79.734 78.738 78.519 26 82.816 82.620 82.347 82.018 80.781 80.006 79.491 78.529 78.268 28 86.751 86.553 86.733 87.321 85.052 85.040 84.449 82.812 83.387 30 94.004 93.806 94.825 97.107 92.961 94.360 93.630 90.745 92.875 32 84.642 84.444 84.380 84.474 82.754 82.330 81.778 80.504 80.628 145 Table C7b. River Hydraulic Heads at Finite Difference Nodes for a 2-Day Time-Step Simulation River Head (ft) Time (days) 9 10 11 12 13 14 15 16 17 0 80.294 82.859 80.659 79.361 79.329 80.275 76.890 77.932 76.528 2 76.192 76.454 75.492 74.866 74.568 74.531 72.976 73.113 72.308 4 100.350 114.174 105.925 101.347 102.619 108.415 96.062 101.620 97.300 6 83.318 87.580 84.469 82.676 82.840 84.515 79.779 81.499 79.654 8 81.142 84.182 81.727 80.291 80.314 81.464 77.700 78.932 77.404 10 82.488 86.284 83.423 81.766 81.877 83.352 78.986 80.520 78.796 12 91.844 100.892 95.209 92.022 92.741 96.478 87.929 91.569 88.485 14 83.567 87.968 84.781 82.948 83.129 84.864 80.016 81.792 79.912 16 78.601 80.215 78.526 77.505 77.363 77.902 75.273 75.938 74.780 18 78.763 80.467 78.730 77.683 77.551 78.128 75.427 76.128 74.947 20 78.464 80.001 78.353 77.355 77.204 77.709 75.142 75.776 74.639 22 77.442 78.405 77.066 76.235 76.018 76.278 74.167 74.576 73.587 24 77.052 77.797 76.575 75.808 75.565 75.733 73.795 74.118 73.187 26 76.860 77.497 76.334 75.598 75.343 75.465 73.612 73.894 72.991 28 80.777 83.613 81.267 79.891 79.890 80.952 77.351 78.501 77.027 30 88.039 94.951 90.416 87.851 88.322 91.139 84.292 87.074 84.543 32 78.666 80.316 78.608 77.576 77.438 77.992 75.334 76.013 74.846 146 Table C7c. River Hydraulic Heads at Finite Difference Nodes for a 2-Day Time-Step Simulation River Head (ft) Time (days) 18 19 20 21 22 23 24 25 26 0 75.558 77.141 75.858 76.196 73.612 76.972 75.633 72.538 70.842 2 71.374 71.481 70.705 70.592 69.301 70.351 69.670 68.298 66.899 4 96.160 105.193 101.403 104.065 95.056 109.926 105.320 93.655 90.480 6 78.659 81.357 79.697 80.380 76.832 81.919 80.090 75.707 73.789 8 76.427 78.322 76.934 77.368 74.514 78.358 76.882 73.426 71.668 10 77.808 80.200 78.643 79.231 75.947 80.561 78.866 74.837 72.980 12 87.417 93.281 90.556 92.227 85.948 95.928 92.709 84.684 82.138 14 78.914 81.704 80.013 80.725 77.097 82.327 80.457 75.969 74.033 16 73.825 74.790 73.717 73.864 71.818 74.216 73.151 70.772 69.200 18 73.990 75.014 73.922 74.086 71.989 74.479 73.387 70.940 69.357 20 73.685 74.600 73.544 73.675 71.673 73.994 72.950 70.629 69.067 22 72.642 73.189 72.260 72.278 70.598 72.343 71.463 69.572 68.084 24 72.246 72.653 71.772 71.749 70.190 71.717 70.900 69.172 67.712 26 72.051 72.391 71.533 71.489 69.991 71.411 70.624 68.975 67.530 28 76.053 77.814 76.471 76.863 74.125 77.761 76.344 73.043 71.312 30 83.507 87.957 85.707 86.936 81.876 89.671 87.073 80.675 78.409 32 73.891 74.879 73.798 73.952 71.885 74.320 73.244 70.838 69.262 147 Table C7d. River Hydraulic Heads at Finite Difference Nodes for a 2-Day Time-Step Simulation River Head (ft) Time (days) 27 28 29 30 31 32 33 34 35 0 72.603 69.907 70.555 67.540 69.021 67.764 66.683 64.530 66.085 2 67.256 65.915 65.617 63.899 63.728 62.495 61.790 60.633 60.969 4 99.260 89.816 95.443 85.894 95.865 94.588 91.647 84.452 92.922 6 76.603 72.894 74.280 70.287 73.037 71.776 70.416 67.509 70.082 8 73.724 70.744 71.598 68.309 70.146 68.888 67.729 65.365 67.204 10 75.505 72.074 73.257 69.532 71.934 70.675 69.391 66.692 68.984 12 87.936 81.358 84.859 78.089 84.449 83.179 81.029 75.979 81.490 14 76.934 73.141 74.588 70.514 73.370 72.108 70.726 67.756 70.415 16 70.376 68.244 68.487 66.015 66.792 65.538 64.611 62.878 63.877 18 70.588 68.402 68.684 66.160 67.004 65.750 64.808 63.034 64.086 20 70.195 68.109 68.320 65.891 66.611 65.358 64.443 62.743 63.699 22 68.862 67.113 67.087 64.982 65.283 64.033 63.210 61.760 62.394 24 68.357 66.736 66.622 64.639 64.782 63.533 62.745 61.389 61.905 26 68.109 66.552 66.395 64.472 64.538 63.289 62.518 61.207 61.666 28 73.241 70.383 71.148 67.977 69.660 68.403 67.277 65.005 66.720 30 82.874 77.578 80.132 74.603 79.349 78.083 76.286 72.194 76.389 32 70.459 68.306 68.564 66.071 66.874 65.620 64.687 62.938 63.958 148 Table C8. Cross Section 1 Aquifer Heads in Layer 1 from 4-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 1 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 110.628 106.773 99.929 88.751 80.288 85.421 91.612 97.677 103.902 4 110.650 106.890 100.671 92.922 95.969 89.857 92.854 98.002 103.985 8 110.677 106.992 100.867 91.573 83.070 88.583 93.153 98.249 104.088 12 110.714 107.120 101.254 92.999 90.663 90.114 93.767 98.550 104.221 16 110.748 107.195 101.207 91.375 81.787 88.444 93.648 98.688 104.331 20 110.770 107.198 100.881 89.583 77.935 86.497 93.070 98.627 104.381 24 110.777 107.139 100.451 88.179 76.024 84.917 92.339 98.419 104.365 28 110.773 107.065 100.189 88.120 78.567 84.763 91.909 98.206 104.312 32 110.771 107.046 100.344 89.917 85.616 86.636 92.190 98.167 104.279 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction. 149 Table C9. Cross Section 1 Aquifer Heads in Layer 2 from 4-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 2 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 101.547 101.651 101.765 101.890 102.039 102.170 102.238 102.256 102.355 4 101.650 101.759 101.877 102.009 102.163 102.292 102.358 102.371 102.462 8 101.642 101.749 101.865 101.993 102.143 102.274 102.342 102.357 102.452 12 101.688 101.797 101.916 102.047 102.200 102.330 102.397 102.410 102.501 16 101.666 101.774 101.890 102.019 102.169 102.300 102.368 102.383 102.477 20 101.626 101.732 101.847 101.973 102.122 102.254 102.322 102.340 102.437 24 101.585 101.690 101.803 101.928 102.075 102.207 102.276 102.295 102.395 28 101.573 101.677 101.791 101.916 102.064 102.196 102.264 102.283 102.383 32 101.604 101.711 101.826 101.954 102.104 102.235 102.303 102.320 102.417 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction. 150 Table C10. Cross Section 2 Aquifer Heads in Layer 1 from 4-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 95.616 90.014 83.018 76.919 72.604 69.927 79.385 87.997 94.995 100.670 4 95.675 90.250 83.953 80.534 87.416 80.988 82.539 88.744 95.159 100.706 8 95.764 90.467 84.273 79.680 75.254 71.902 81.736 88.986 95.309 100.758 12 95.885 90.740 84.827 81.018 82.405 77.249 82.881 89.417 95.496 100.830 16 96.002 90.916 84.856 79.797 74.037 70.996 81.779 89.424 95.615 100.899 20 96.083 90.943 84.505 78.259 70.409 68.289 80.432 89.134 95.631 100.948 24 96.110 90.839 83.980 76.937 68.635 66.955 79.299 88.716 95.559 100.965 28 96.098 90.699 83.623 76.708 70.990 68.714 79.132 88.439 95.463 100.961 32 96.086 90.661 83.778 78.143 77.630 73.681 80.411 88.575 95.441 100.957 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 151 Table C11. Cross Section 2 Aquifer Heads in Layer 2 from 4-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 2 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 98.268 98.074 97.901 97.677 97.583 97.596 97.676 97.789 97.873 97.978 4 98.368 98.180 98.016 97.799 97.703 97.694 97.750 97.840 97.914 98.017 8 98.348 98.155 97.984 97.758 97.655 97.657 97.726 97.827 97.905 98.009 12 98.392 98.202 98.034 97.811 97.707 97.700 97.758 97.850 97.923 98.026 16 98.366 98.172 97.999 97.769 97.663 97.665 97.734 97.835 97.912 98.016 20 98.329 98.133 97.957 97.727 97.623 97.632 97.709 97.817 97.898 98.003 24 98.294 98.097 97.921 97.691 97.591 97.605 97.687 97.802 97.886 97.991 28 98.287 98.091 97.916 97.689 97.592 97.605 97.686 97.801 97.884 97.989 32 98.321 98.129 97.958 97.735 97.638 97.642 97.713 97.818 97.898 98.002 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 152 Table C12. Leakage Rates from 4-Day Time-Step Simulation Leakage Rate (ft3/d) Time (days) Reach14 Reach 27 Reach 28 Total 0 -375.201 -174.686 -821.384 -21216.420 4 1029.300 2587.630 1069.640 299930.497 8 -964.381 -1627.310 -2014.150 -250130.969 12 351.587 1215.300 50.630 123104.650 16 -836.333 -1273.370 -1739.830 -185140.890 20 -726.143 -886.831 -1379.340 -109646.680 24 -672.762 -715.988 -1201.160 -69520.910 28 -319.623 -5.775 -670.439 22445.224 32 195.845 959.140 -22.168 123464.396 153 Table C13a. River Discharge at Finite Difference Nodes for a 4-Day Time-Step Simulation River Discharge (m3/s) Time (days) 0 1 2 3 4 5 6 7 8 0 77.300 77.348 77.416 77.456 77.982 77.991 78.008 78.062 78.181 4 412.013 412.062 412.130 412.170 412.696 412.705 412.722 412.776 412.895 8 120.770 120.810 120.871 120.894 121.409 121.400 121.398 121.451 121.568 12 275.662 275.717 275.790 275.843 276.378 276.400 276.431 276.486 276.606 16 99.745 99.790 99.854 99.887 100.409 100.410 100.419 100.472 100.590 20 46.722 46.775 46.847 46.896 47.429 47.447 47.474 47.529 47.649 24 26.603 26.654 26.724 26.769 27.299 27.312 27.335 27.389 27.509 28 54.417 54.467 54.536 54.578 55.107 55.118 55.138 55.192 55.311 32 167.280 167.327 167.394 167.432 167.957 167.963 167.977 168.031 168.149 154 Table C13b. River Discharge at Finite Difference Nodes for a 4-Day Time-Step Simulation River Discharge (m3/s) Time (days) 9 10 11 12 13 14 15 16 17 0 78.192 78.206 78.220 78.226 78.240 78.393 78.396 78.761 78.908 4 412.907 412.921 412.935 412.942 412.956 413.110 413.112 413.478 413.625 8 121.578 121.592 121.606 121.610 121.624 121.777 121.779 122.145 122.292 12 276.619 276.633 276.648 276.656 276.671 276.824 276.827 277.193 277.340 16 100.601 100.615 100.629 100.635 100.649 100.802 100.805 101.170 101.317 20 47.661 47.675 47.690 47.698 47.712 47.865 47.868 48.234 48.382 24 27.521 27.535 27.550 27.557 27.571 27.725 27.728 28.093 28.241 28 55.323 55.337 55.352 55.359 55.373 55.527 55.530 55.895 56.043 32 168.161 168.175 168.189 168.196 168.210 168.363 168.366 168.732 168.879 155 Table C13c. River Discharge at Finite Difference Nodes for a 4-Day Time-Step Simulation River Discharge (m3/s) Time (days) 18 19 20 21 22 23 24 25 26 0 78.945 79.611 79.642 79.982 80.010 80.069 80.089 80.134 80.137 4 413.662 414.328 414.359 414.699 414.728 414.787 414.807 414.853 414.856 8 122.328 122.993 123.024 123.363 123.392 123.450 123.470 123.515 123.517 12 277.377 278.043 278.075 278.415 278.444 278.504 278.524 278.571 278.574 16 101.354 102.019 102.050 102.390 102.418 102.477 102.497 102.543 102.545 20 48.419 49.084 49.116 49.456 49.485 49.545 49.565 49.612 49.615 24 28.278 28.943 28.975 29.315 29.344 29.403 29.423 29.470 29.473 28 56.080 56.745 56.777 57.117 57.145 57.205 57.225 57.271 57.274 32 168.916 169.582 169.613 169.953 169.981 170.041 170.061 170.107 170.109 156 Table C13d. River Discharge at Finite Difference Nodes for a 4-Day Time-Step Simulation River Discharge (m3/s) Time (days) 27 28 29 30 31 32 33 34 35 0 80.211 80.242 81.366 81.386 81.485 81.615 81.632 81.647 83.479 4 414.930 414.961 416.086 416.106 416.206 416.336 416.354 416.368 418.201 8 123.590 123.621 124.744 124.764 124.859 124.989 125.006 125.019 126.851 12 278.648 278.680 279.806 279.826 279.929 280.059 280.078 280.093 281.925 16 102.619 102.650 103.774 103.794 103.891 104.022 104.039 104.053 105.885 20 49.689 49.721 50.846 50.866 50.968 51.099 51.117 51.132 52.964 24 29.547 29.578 30.703 30.723 30.825 30.955 30.973 30.988 32.820 28 57.348 57.380 58.505 58.525 58.626 58.756 58.774 58.789 60.621 32 170.183 170.215 171.339 171.359 171.459 171.589 171.607 171.621 173.453 157 Table C14a. River Hydraulic Head at Finite Difference Nodes for a 4-Day Time-Step Simulation River Head (ft) Time (days) 0 1 2 3 4 5 6 7 8 0 86.262 86.070 86.194 86.670 84.526 84.420 83.838 82.284 82.756 4 97.451 97.253 98.671 101.759 96.724 98.795 98.000 94.522 97.393 8 88.253 88.055 88.408 89.346 86.687 86.965 86.345 84.449 85.345 12 93.669 93.471 94.452 96.656 92.598 93.932 93.210 90.382 92.441 16 87.336 87.139 87.386 88.110 85.689 85.790 85.188 83.450 84.150 20 84.581 84.385 84.314 84.398 82.694 82.262 81.713 80.448 80.561 24 83.210 83.014 82.786 82.550 81.207 80.509 79.987 78.957 78.779 28 85.039 84.842 84.824 85.014 83.190 82.846 82.288 80.945 81.155 32 90.078 89.880 90.444 91.809 88.678 89.313 88.658 86.449 87.736 158 Table C14b. River Hydraulic Head at Finite Difference Nodes for a 4-Day Time-Step Simulation River Head (ft) Time (days) 9 10 11 12 13 14 15 16 17 0 80.294 82.859 80.659 79.361 79.329 80.275 76.890 77.932 76.528 4 91.496 100.350 94.771 91.641 92.337 95.990 87.597 91.158 88.125 8 82.276 85.953 83.155 81.533 81.630 83.054 78.783 80.269 78.576 12 87.707 94.432 89.997 87.487 87.937 90.673 83.974 86.681 84.199 16 81.361 84.524 82.003 80.531 80.568 81.771 77.909 79.190 77.630 20 78.615 80.236 78.544 77.521 77.379 77.921 75.286 75.954 74.794 24 77.251 78.108 76.827 76.027 75.797 76.012 73.985 74.352 73.392 28 79.069 80.946 79.116 78.019 77.907 78.558 75.720 76.489 75.263 32 84.106 88.810 85.461 83.540 83.755 85.621 80.532 82.429 80.470 159 Table C14c. River Hydraulic Head at Finite Difference Nodes for a 4-Day Time-Step Simulation River Head (ft) Time (days) 18 19 20 21 22 23 24 25 26 0 75.558 77.141 75.858 76.196 73.612 76.972 75.633 72.538 70.842 4 87.060 92.795 90.113 91.744 85.576 95.356 92.194 84.318 81.797 8 77.589 79.903 78.373 78.937 75.721 80.212 78.552 74.613 72.772 12 83.166 87.492 85.284 86.475 81.521 89.126 86.582 80.325 78.084 16 76.651 78.627 77.211 77.670 74.746 78.715 77.203 73.654 71.880 20 73.839 74.809 73.735 73.883 71.832 74.239 73.171 70.786 69.213 24 72.449 72.927 72.022 72.020 70.399 72.037 71.188 69.376 67.902 28 74.304 75.439 74.308 74.507 72.312 74.976 73.835 71.259 69.653 32 79.468 82.458 80.700 81.474 77.673 83.213 81.255 76.536 74.560 160 Table C14d. River Hydraulic Head at Finite Difference Nodes for a 4-Day Time-Step Simulation River Head (ft) Time (days) 27 28 29 30 31 32 33 34 35 0 72.603 69.907 70.555 67.540 69.021 67.764 66.683 64.530 66.085 4 87.473 81.012 84.428 77.770 83.983 82.713 80.596 75.633 81.024 8 75.223 71.864 72.994 69.338 71.650 70.391 69.127 66.481 68.701 12 82.433 77.248 79.720 74.299 78.905 77.640 75.874 71.865 75.945 16 74.013 70.959 71.866 68.507 70.434 69.176 67.997 65.579 67.490 20 70.394 68.257 68.504 66.027 66.809 65.556 64.627 62.890 63.894 24 68.615 66.929 66.860 64.815 65.038 63.788 62.982 61.578 62.153 28 70.990 68.702 69.056 66.434 67.405 66.151 65.181 63.332 64.483 32 77.650 73.676 75.256 71.007 74.090 72.829 71.396 68.291 71.134 161 Table C15. Cross Section 1 Aquifer Heads in Layer 1 from 8-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 1 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 110.628 106.773 99.929 88.751 80.288 85.421 91.612 97.677 103.902 8 110.679 106.983 100.840 92.322 90.265 89.304 93.106 98.214 104.088 16 110.748 107.180 101.249 92.312 86.656 89.400 93.727 98.649 104.313 24 110.770 107.134 100.601 88.943 77.037 85.758 92.607 98.437 104.339 32 110.780 107.106 100.530 89.714 82.474 86.501 92.491 98.335 104.328 Table C16. Cross Section 1 Aquifer Heads in Layer 2 from 8-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 2 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 101.547 101.651 101.765 101.890 102.039 102.170 102.238 102.256 101.547 8 101.652 101.760 101.878 102.008 102.160 102.291 102.357 102.372 101.652 16 101.681 101.790 101.907 102.037 102.189 102.320 102.387 102.401 101.681 24 101.603 101.708 101.822 101.948 102.096 102.228 102.296 102.314 101.603 32 101.613 101.719 101.834 101.962 102.112 102.243 102.311 102.328 101.613 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction . 162 Table C17. Cross Section 2 Aquifer Heads in Layer 1 from 8-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 95.616 90.014 83.018 76.919 72.604 69.927 79.385 87.997 94.995 100.670 8 95.773 90.452 84.234 80.216 82.030 76.973 82.224 88.956 95.298 100.763 16 95.998 90.886 84.889 80.522 78.633 74.429 82.429 89.456 95.595 100.898 24 96.082 90.826 84.141 77.598 69.571 67.660 79.879 88.852 95.553 100.951 32 96.116 90.784 84.040 78.119 74.665 71.462 80.362 88.779 95.527 100.976 Table C18. Cross Section 2 Aquifer Heads in Layer 2 from 8-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 2 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 98.268 98.074 97.901 97.677 97.583 97.596 97.676 97.789 97.873 97.978 8 98.364 98.174 98.006 97.784 97.685 97.681 97.742 97.837 97.912 98.016 16 98.384 98.192 98.022 97.795 97.690 97.686 97.749 97.844 97.919 98.023 24 98.310 98.113 97.938 97.708 97.606 97.618 97.697 97.809 97.891 97.996 32 98.324 98.130 97.957 97.731 97.631 97.638 97.711 97.818 97.898 98.003 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 163 Table C19. Leakage Rates from 8-Day Time-Step Simulation Leakage Rate (ft3/d) Time (days) Reach14 Reach 27 Reach 28 Total 0 -375.201 -174.686 -821.384 -21216.420 8 239.400 941.619 -129.422 89179.734 16 -279.300 -135.934 -930.565 -48144.746 24 -789.873 -1019.400 -1455.530 -123170.194 32 -172.230 221.782 -563.729 31542.606 164 Table C20a. River Discharge at Finite Difference Nodes for 8-Day Time-Step Simulation River Discharge (m3/s) Time (days) 0 1 2 3 4 5 6 7 8 0 77.300 77.348 77.416 77.456 77.982 77.991 78.008 78.062 78.181 8 266.392 266.440 266.508 266.548 267.075 267.083 267.101 267.155 267.274 16 187.704 187.749 187.815 187.849 188.372 188.375 188.385 188.439 188.557 24 36.663 36.712 36.781 36.822 37.350 37.361 37.380 37.434 37.553 32 110.849 110.900 110.970 111.015 111.546 111.560 111.584 111.638 111.758 Table C20b. River Discharge at Finite Difference Nodes for 8-Day Time-Step Simulation River Discharge (m3/s) Time (days) 9 10 11 12 13 14 15 16 17 0 78.192 78.206 78.220 78.226 78.240 78.393 78.396 78.761 78.908 8 267.285 267.299 267.314 267.321 267.335 267.488 267.491 267.856 268.004 16 188.568 188.582 188.597 188.602 188.617 188.769 188.772 189.138 189.285 24 37.565 37.579 37.593 37.600 37.615 37.768 37.771 38.136 38.283 32 111.770 111.784 111.798 111.806 111.820 111.974 111.977 112.342 112.490 165 Table C20c. River Discharge at Finite Difference Nodes for 8-Day Time-Step Simulation River Discharge (m3/s) Time (days) 18 19 20 21 22 23 24 25 26 0 78.945 79.611 79.642 79.982 80.010 80.069 80.089 80.134 80.137 8 268.041 268.706 268.738 269.078 269.106 269.166 269.186 269.232 269.234 16 189.322 189.987 190.018 190.358 190.386 190.446 190.465 190.511 190.514 24 38.320 38.986 39.017 39.357 39.386 39.445 39.465 39.511 39.514 32 112.527 113.192 113.224 113.564 113.593 113.653 113.673 113.719 113.722 Table C20d. River Discharge at Finite Difference Nodes for 8-Day Time-Step Simulation River Discharge (m3/s) Time (days) 27 28 29 30 31 32 33 34 35 0 80.211 80.242 81.366 81.386 81.485 81.615 81.632 81.647 83.479 8 269.308 269.340 270.464 270.484 270.585 270.715 270.732 270.747 272.579 16 190.587 190.618 191.743 191.762 191.861 191.991 192.008 192.022 193.854 24 39.588 39.620 40.744 40.764 40.865 40.995 41.012 41.027 42.859 32 113.796 113.828 114.953 114.973 115.075 115.205 115.223 115.238 117.070 166 Table C21a. River Hydraulic Head at Finite Difference Nodes for 8-Day Time-Step Simulation River Head (ft) Time (days) 0 1 2 3 4 5 6 7 8 0 86.262 86.070 86.194 86.670 84.526 84.420 83.838 82.284 82.756 8 93.388 93.189 94.137 96.275 92.289 93.568 92.851 90.072 92.070 16 90.813 90.615 91.264 92.801 89.480 90.257 89.589 87.253 88.698 24 83.934 83.737 83.592 83.524 81.990 81.432 80.895 79.741 79.716 32 87.829 87.631 87.936 88.776 86.227 86.425 85.814 83.991 84.797 Table C21b. River Hydraulic Head at Finite Difference Nodes for 8-Day Time-Step Simulation River Head (ft) Time (days) 9 10 11 12 13 14 15 16 17 0 80.294 82.859 80.659 79.361 79.329 80.275 76.890 77.932 76.528 8 87.422 93.988 89.639 87.175 87.606 90.274 83.702 86.346 83.904 16 84.842 89.959 86.388 84.346 84.610 86.653 81.235 83.298 81.232 24 77.968 79.227 77.730 76.812 76.629 77.016 74.669 75.194 74.129 32 81.857 85.298 82.627 81.074 81.143 82.466 78.383 79.775 78.143 167 Table C21c. River Hydraulic Head at Finite Difference Nodes for 8-Day Time-Step Simulation River Head (ft) Time (days) 18 19 20 21 22 23 24 25 26 0 75.558 77.141 75.858 76.196 73.612 76.972 75.633 72.538 70.842 8 82.874 87.094 84.922 86.080 81.217 88.658 86.160 80.026 77.806 16 80.223 83.486 81.636 82.495 78.459 84.420 82.343 77.310 75.280 24 73.179 73.915 72.921 72.997 71.151 73.192 72.228 70.116 68.590 32 77.160 79.319 77.841 78.357 75.274 79.527 77.934 74.174 72.364 Table C21d. River Hydraulic Head at Finite Difference Nodes for 8-Day Time-Step Simulation River Head (ft) Time (days) 27 28 29 30 31 32 33 34 35 0 72.603 69.907 70.555 67.540 69.021 67.764 66.683 64.530 66.085 8 82.055 76.966 79.367 74.038 78.524 77.259 75.519 71.582 75.564 16 78.627 74.405 76.168 71.679 75.073 73.810 72.310 69.021 72.115 24 69.547 67.625 67.720 65.449 65.965 64.713 63.843 62.264 63.062 32 74.669 71.450 72.478 68.958 71.094 69.835 68.610 66.068 68.147 168 Table C22. Cross Section 1 Aquifer Heads in Layer 1 from 16-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 1 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 110.628 106.773 99.929 88.751 80.288 85.421 91.612 97.677 103.902 16 110.747 107.148 101.139 92.377 88.530 89.423 93.546 98.555 104.290 32 110.777 107.112 100.580 89.517 80.031 86.332 92.579 98.366 104.322 Table C23. Cross Section 1 Aquifer Heads in Layer 2 from 16-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 1, Layer 2 Time (days) -4 -3 -2 -1* Reach 14 1 2 3 4 0 101.547 101.651 101.765 101.89 102.039 102.17 102.238 102.256 102.355 16 101.677 101.785 101.903 102.033 102.185 102.316 102.383 102.397 102.489 32 101.612 101.718 101.832 101.959 102.108 102.239 102.307 102.325 102.423 *Negative values indicate the distance in MODFLOW cells from Reach 14 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 14 in the northeastern direction. 169 Table C24. Cross Section 2 Aquifer Heads in Layer 1 from 16-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 95.616 90.014 83.018 76.919 72.604 69.927 79.385 87.997 94.995 100.670 16 95.996 90.826 84.738 80.492 80.391 75.740 82.418 89.339 95.555 100.898 32 96.105 90.793 84.101 78.005 72.381 69.753 80.251 88.828 95.536 100.972 Table C25. Cross Section 2 Aquifer Heads in Layer 1 from 16-Day Time-Step Simulation Aquifer Heads (ft) in Cross Section 2, Layer 1 Time (days) -4 -3 -2 -1* Reach 27 Reach 28 1 2 3 4 0 98.268 98.074 97.901 97.677 97.583 97.596 97.676 97.789 97.873 97.978 16 98.382 98.191 98.022 97.797 97.693 97.689 97.750 97.845 97.919 98.023 32 98.321 98.126 97.952 97.724 97.623 97.631 97.707 97.816 97.897 98.002 *Negative values indicate the distance in MODFLOW cells from Reach 27 in the southwestern direction. Positive values indicate the distance in MODFLOW cells from Reach 28 in the northeastern direction. 170 Table C26. Leakage Rates from 16-Day Time-Step Simulation Leakage Rate (ft3/d) Time (days) Reach14 Reach 27 Reach 28 Total 0 -375.201 -174.686 -821.384 -21216.420 16 2.850 453.130 -495.138 28816.784 32 -525.836 -509.674 -1100.140 -64873.732 171 Table C27a. River Discharge at Finite Difference Nodes for 16-Day Time-Step Simulation River Discharge (m3/s) Time (days) 0 1 2 3 4 5 6 7 8 0 77.300 77.348 77.416 77.456 77.982 77.991 78.008 78.062 78.181 16 227.048 227.096 227.164 227.204 227.731 227.739 227.756 227.810 227.929 32 73.756 73.803 73.870 73.907 74.433 74.439 74.453 74.507 74.626 Table C27b. River Discharge at Finite Difference Nodes for 16-Day Time-Step Simulation River Discharge (m3/s) Time (days) 9 10 11 12 13 14 15 16 17 0 78.192 78.206 78.220 78.226 78.240 78.393 78.396 78.761 78.908 16 227.941 227.955 227.970 227.977 227.991 228.144 228.147 228.512 228.660 32 74.637 74.651 74.665 74.672 74.686 74.839 74.842 75.207 75.354 172 Table C27c. River Discharge at Finite Difference Nodes for 16-Day Time-Step Simulation River Discharge (m3/s) Time (days) 18 19 20 21 22 23 24 25 26 0 78.945 79.611 79.642 79.982 80.010 80.069 80.089 80.134 80.137 16 228.697 229.362 229.393 229.733 229.762 229.822 229.842 229.888 229.890 32 75.391 76.057 76.088 76.428 76.456 76.516 76.536 76.581 76.584 Table C27d. River Discharge at Finite Difference Nodes for 16-Day Time-Step Simulation River Discharge (m3/s) Time (days) 27 28 29 30 31 32 33 34 35 0 80.211 80.242 81.366 81.386 81.485 81.615 81.632 81.647 83.479 16 229.964 229.996 231.120 231.140 231.241 231.371 231.388 231.403 233.235 32 76.658 76.689 77.814 77.833 77.933 78.063 78.080 78.095 79.927 173 Table C28a. River Hydraulic Head at Finite Difference Nodes for 16-Day Time-Step Simulation River Head (ft) Time (days) 0 1 2 3 4 5 6 7 8 0 86.262 86.070 86.194 86.670 84.526 84.420 83.838 82.284 82.756 16 92.145 91.947 92.751 94.599 90.933 91.971 91.277 88.712 90.443 32 86.088 85.890 85.993 86.427 84.330 84.189 83.611 82.087 82.520 Table C28b. River Hydraulic Head at Finite Difference Nodes for 16-Day Time-Step Simulation River Head (ft) Time (days) 9 10 11 12 13 14 15 16 17 0 80.294 82.859 80.659 79.361 79.329 80.275 76.890 77.932 76.528 16 86.177 92.045 88.071 85.811 86.161 88.527 82.512 84.875 82.615 32 80.114 82.578 80.432 79.164 79.120 80.023 76.718 77.720 76.342 174 Table C28c. River Hydraulic Head at Finite Difference Nodes for 16-Day Time-Step Simulation River Head (ft) Time (days) 18 19 20 21 22 23 24 25 26 0 75.558 77.141 75.858 76.196 73.612 76.972 75.633 72.538 70.842 16 81.595 85.353 83.336 84.350 79.886 86.613 84.318 78.715 76.587 32 75.373 76.891 75.630 75.947 73.421 76.679 75.369 72.349 70.667 Table C28d. River Hydraulic Head at Finite Difference Nodes for 16-Day Time-Step Simulation River Head (ft) Time (days) 27 28 29 30 31 32 33 34 35 0 72.603 69.907 70.555 67.540 69.021 67.764 66.683 64.530 66.085 16 80.400 75.730 77.823 72.899 76.859 75.594 73.970 70.346 73.899 32 72.366 69.730 70.334 67.377 68.783 67.526 66.461 64.354 65.848 175 References Abbott, M. B., J. C. Bathurst, J.A. Cunge, P.E. O'Connell and J. Rasmussen. ?An Introduction to the European Hydrological System- Systeme Hydrologique Europeen, 'SHE', 1: History and Philosophy of a physically- based , distributed modeling system.? Journal of Hydrology. 87 (1986): 45-59. Bathurst, J.C. "Sensitivity Analysis of the Systeme Hydrologique Europeen for an Updated Catchment." Journal of Hydrology. 87 (1986): 103-123. Charbeneau, R. Groundwater Hydraulics and Pollutant Transport. Upper Saddle River: Prentice Hall, 2000. Chow, V. T., D. R. Maidment, and L. W. Mays. Applied Hydrology. New York: McGraw Hill, 1988. Danish Hydraulic Institute. "MIKE-11 - A Modeling System for Rivers and Channels" Online. Available: http://www.dhi.dk/mike11/index.htm. Accessed April 2, 2000. Danish Hydraulic Institute. "MIKE-SHE - An Integrated Hydrological Modeling System." Online. Available: http://www.dhi.dk/mikeshe/index.htm. Accessed April 2, 2000. Dutton, A. R. Groundwater Availability in the Carrizo-Wilcox Aquifer in Central Texas- Numerical Simulations of 2000 through 2050 Withdrawal Projections. Report No. 256. Austin, Tex: Bureau of Economic Geology, 1999. Hinaman, K. C. ?Use of a Geographic Information System to Assemble Input- Data Sets for a Finite Difference Model of Groundwater Flow? Water Resource Bulletin, American Water Resources Association. 29 (1993) :401-406. Hubert, M. and B. Bullock. ?Senate Bill 1, The First Big and Bold Step Toward Meeting Texas?s Future Water Needs? 30, Texas Technical Law Review. 53 (1999). HydroGeologic, Inc. "MODFLOW-SURFACT Volume III: Surface Water Flow Modules." Herndon, VA. (Software Documentation) 176 Lighthill, M. J., and G.B. Whitham. "On Kinematic Waves, I: Flood Movement in Long Rivers." Proc. R. Soc. London A. 229 (1955):281-316. Lusk, Stephanie E. Hayes ?Texas Groundwater: Reconciling the Rule of Capture with the Environmental and Community Demands,? 30 St. Mary?s Law Journal 305, 1998. Mace, R. E., and Mullican, W. F., III, 2000, The past, present, and future of groundwater availability modeling in Texas: Southwest Focus Ground Water Conference, National Ground Water Association, p. 35-36. Mueller, F. A. and J. W. Male. ?A Management Model for Specification of Groundwater Withdrawal Permits.? Water Resources Research. 29 (1993): 1359-1368. Olivera, Francisco. "Calculation of Hydrologic Parameters using CRWR PrePro." Presented in October 1998 at the CE394-K3 GIS in Water Resources class at the University of Texas at Austin. Austin, Texas. Online. Available: http://www.ce.utexas.edu/prof/olivera/prepro/prepro.htm. Accessed October 13, 1999. Orzol, L. L. and T. S. McGrath. ?Summary of Modifications of the U.S. Geological Survey Modular Finite Difference, Groundwater Flow Model to Read and Write Geographic Information System Files.? Water Resources Bulletin, American Water Resources Association. 29 (1993):843-846. Perkins, S. P. and M. Sophocleous. "Development of a Comprehensive Watershed Model Applied to Study Stream Yield under Drought Conditions." Ground Water. 37 (1999): 418-426. Richards, C. J., h. Roaza and R. M. Roaza. "Integrating Geographic Information system and MODFLOW for Ground Water Resources Assessments." Water Resources Bulletin, American Water Resources Association. 29 (1993): 847-853. Refsgaard, J. C. ?Parameterization, Calibration and Validation of Distributed Hydrological Models.? Journal of Hydrology. 198 (1997): 69-97. Refsgaard, J. C. and J. Knudsen ?Operational validation and intercomparison of different types of hydrological models? Water Resources Research. 32 (1996):2189-2202. 177 Romanek, Andrew "Surface Water Modeling at Marcus Hook." Presented on April 23, 1998. Online. Available: http://www.ce.utexas.edu/prof/maidment/grad/romanek/gisproject/gis/sld0 01.htm. Accessed March 17, 2000. Sophocleous, M. ?Evaluation of Simplified Stream Aquifer Depletion Models for Water Rights Administration.? Ground Water. 33 (1995): 579-588. Texas Water Code 16.012, 75th Legislature (1997). Online. Available: http://capitol.tlc.state.tx.us/statutes/codes/WA000010.html. Accessed April 14, 2000. Texas Water Code 36.311, 75th Legislature (1997). Online. Available: http://capitol.tlc.state.tx.us/statutes/codes/WA000025.html. Accessed April 14, 2000. Vance, Benjamin R. ?Total Aquifer Management: A New Approach to Groundwater Protection.? 30, University of San Francisco Law Review. 803 (1996). Vieira, J. H. ?Conditions governing the use of the approximations for the Saint- Venant Equations for Shallow Surface Water Flow.? Journal of Hydrology. 60 (1983): 43-58. 178 Vita Shiva Niazi was born in Mashad, Iran on September 21, 1975. She moved to Berkeley, CA in August of 1979 following the Revolution in Iran and has resided in America ever since. She is the youngest of four children. Her parents, Mansour and Fereshteh Niazi, sister Jaleh, and brother Kaveh still live in the San Francisco Bay Area. Her third sibling, Kayvan is currently living in Los Angeles, CA. Shiva attended Berkeley High School from 1989 through 1993. She went on to receive her B.S. in Environmental Engineering Science in 1997 from the University of California at Berkeley. Following her undergraduate work, she worked for an environmental consultant, SOMA Engineering, Inc. for six months before beginning her graduate work at the University of Texas at Austin in August 1998. Permanent address: 925 Hilldale Ave. Berkeley, CA 94708 This thesis was typed by the author.