Therefore, in Nebraska, groundwater is the most important natural resource and needs to .... Ogallala Group (Nebraska Department of Natural Resources, 2005). ...... Mead. Yutan. MEADAGROFARM. 6804900. Nebraska. Legend. City ...... Changes in hyporheic exchange flow following experimental wood removal in a small ...
MODELING SURFACE WATER AND GROUNDWATER INTERACTIONS IN AGRICULTURAL AREAS
by
Gengxin Ou
A DISSERTATION
Presented to the Faculty of The Graduate College at the University of Nebraska In Partial Fulfilment of Requirements For the Degree of Doctor of Philosophy
Major: Civil Engineering
Under the Supervision of Ayse Kilic and Xun-Hong Chen
Lincoln, Nebraska May, 2015
ProQuest Number: 3718063
All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion.
ProQuest 3718063 Published by ProQuest LLC (2015). Copyright of the Dissertation is held by the Author. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code Microform Edition © ProQuest LLC. ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106 - 1346
MODELING SURFACE WATER AND GROUNDWATER INTERACTIONS IN AGRICULTURAL AREAS Gengxin Ou, Ph.D. University of Nebraska, 2015 Adviser: Ayse Kilic and Xun-Hong Chen In this dissertation research, three modeling tools have been developed to help quantify the complex processes of integrated surface- and ground- water systems in an area under significant agricultural development of Nebraska. To evaluate irrigation effect on stream depletion for the Lower Platte River (LPR) basin, a regional three-dimensional, transient groundwater-flow model is established with MODFLOW. A simplified and more efficient solution has been developed to estimate the stream depletion rate by re-using the flow coefficients of the baseline run. The flow equation is linearized because the head coefficients become constant between solver iterations. The tool has been successfully applied to the LPR model. The results show that the stream depletion analysis tool can reduce the numerical errors produced in the conventional method and improve the computational efficiency. The second modeling tool is a cross-section based streamflow routing (CSR) package for MODFLOW to simulate the streamflow and the interaction between streams and aquifers for the stream with a width larger than the MODFLOW grid size. In the CSR package, streams are divided into stream segments which are formed by two consecutive cross-sections. A rapid algorithm is used to compute the submerged area of the MODFLOW grid. Streamaquifer seepage is treated as lateral flow in the streamflow routing computation with the Muskingum-Cunge method or mass conservation method. A hypothetical
problem was established to test the capabilities of the CSR package with steadyand transient-state models. The third part of this dissertation aims to quantify the impacts of natural processes and human activities on ground-water dynamics in highly agricultural areas by developing an integrated surface-ground water model (ISGWM). In ISGWM, SWAT and MODFLOW are linked by a soil water module (SWM), which is developed based on a non-iterative solution of the 1D Richards equation. SWM explicitly represents infiltration, soil evaporation, unsaturated water flow, root water update, and lateral drainage. Taking advantage of the simulation capacities of these SWAT, MODFLOW and SWM, ISGWM can simulate the physical hydrologic processes in three domains and their interactions. The model has been successfully applied to the Johnson Creek watershed, which is located within the Lower Platte Basin.
iv DEDICATION I dedicate this dissertation to my family as they have always supported me no matter what situation I was facing.
v ACKNOWLEDGMENTS This dissertation could not have been finished without the help and support from many professors, colleagues, friends and my family. First of all, I would like to thank Drs. Ayse Kilic and Xun-Hong Chen. I especially owe Dr. Ayse Kilic a huge debt of gratitude for her help on my PhD program. She always made time to discuss my progress and give me helpful and quick feedback. I thank Dr. Xun-Hong Chen, for his inspiration, helpful instructions, and tremendous efforts and dedication to guide me through my study and research over the years. I also would like to thank my supervisory committee members, Drs. Shannon Bartelt-Hunt, Yusong Li, and Ashok Samal for their great and enthusiastic help and useful advices on both of my research and the writing of the dissertation. I would also like to thank my colleagues, Dr. Cheng Cheng, Dr. Ruopu Li, Mahesh Pun, Can Liu, Zhaowei Wang. I received help from them time from time. Their friendships and suggestions provide comfort even I was tens of thousands miles from my home. Last but not least, I would like to thank my parents Haiqi Ou, Mingyu Chen and wife Jiong Wu for their endless deep love and good care for me. Their constant support and encouragement helped me successfully complete this challenging venture. My deepest appreciation is expressed to my wife who is always being there for me. Her company gave me strength to tackle all kinds of difficulties.
vi
Contents Contents
vii
List of Figures
xii
List of Tables
xv
1 INTRODUCTION
1
1.1
Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Research summary and objectives . . . . . . . . . . . . . . . . . . . .
4
1.3
Agenda of the dissertation . . . . . . . . . . . . . . . . . . . . . . . . .
7
2 STREAM DEPLETION ANALYSIS IN THE LOWER PLATTE RIVER 8
BASIN 2.1
2.2
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.1.1
Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.1.2
Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.1.3
Previous studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Conceptual groundwater-flow model . . . . . . . . . . . . . . . . . . 16 2.2.1
Hydrostratigraphic units . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1.1
Hydrogeologic settings . . . . . . . . . . . . . . . . . 16
vii 2.2.1.2
2.3
2.2.2
Groundwater levels . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2.3
Groundwater pumping . . . . . . . . . . . . . . . . . . . . . . 26
2.2.4
Recharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.5
Evapotranspiration . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.6
Surface water . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6
Stream network . . . . . . . . . . . . . . . . . . . . . . 29
2.2.6.2
Stream stages . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1
Space and time discretization . . . . . . . . . . . . . . . . . . . 33
2.3.2
Initial condition . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.3
Hydraulic properties . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.4
Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3.4.1
Lateral boundary conditions . . . . . . . . . . . . . . 41
2.3.4.2
Groundwater recharge . . . . . . . . . . . . . . . . . . 41
2.3.4.3
Groundwater evapotranspiration . . . . . . . . . . . 43
2.3.4.4
Streambed seepage . . . . . . . . . . . . . . . . . . . . 45
2.3.4.5
Pumping wells . . . . . . . . . . . . . . . . . . . . . . 46
MODFLOW Packages . . . . . . . . . . . . . . . . . . . . . . . 46
Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.4.1
2.5
2.2.6.1
Development of the numerical model . . . . . . . . . . . . . . . . . . 33
2.3.5 2.4
Elevation calculation . . . . . . . . . . . . . . . . . . . 19
Parameters in model calibration . . . . . . . . . . . . . . . . . 49
Model results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.5.1
Regional water budget
. . . . . . . . . . . . . . . . . . . . . . 50
2.5.2
Groundwater levels . . . . . . . . . . . . . . . . . . . . . . . . . 51
Streamflow depletion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.6.1
Explanation of MODFLOW-SDA . . . . . . . . . . . . . . . . . 55
viii 2.6.2
2.6.3 2.7
Application of MODFLOW-SDA . . . . . . . . . . . . . . . . . 59 2.6.2.1
Hypothetical case . . . . . . . . . . . . . . . . . . . . 60
2.6.2.2
LPR model case . . . . . . . . . . . . . . . . . . . . . . 61
Discussion on MODFLOW-SDA . . . . . . . . . . . . . . . . . 66
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3 A CROSS-SECTION BASED STREAMFLOW ROUTING PACKAGE FOR MODFLOW
69
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.2
Literature review on River and Stream Packages of MODFLOW . . . 71
3.3
Development of the CSR package . . . . . . . . . . . . . . . . . . . . . 77
3.4
3.5
3.6
3.3.1
Overview of the model . . . . . . . . . . . . . . . . . . . . . . . 77
3.3.2
Conceptualization of stream components . . . . . . . . . . . . 78
3.3.3
Computation of stream depth . . . . . . . . . . . . . . . . . . . 79
3.3.4
Computation of intersections of the submerged area . . . . . 83
3.3.5
Computation of the stream-aquifer seepage . . . . . . . . . . . 85
3.3.6
Streamflow routing . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.3.7
Implementation of CSR with MODFLOW . . . . . . . . . . . . 89
Application of the CSR package . . . . . . . . . . . . . . . . . . . . . . 90 3.4.1
Test case: steady-state model . . . . . . . . . . . . . . . . . . . 92
3.4.2
Test case: transient-state model . . . . . . . . . . . . . . . . . . 96
Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 3.5.1
Advantages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.5.2
Assumptions and limitations . . . . . . . . . . . . . . . . . . . 99
3.5.3
Potential future improvements . . . . . . . . . . . . . . . . . . 101
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
ix 4 AN INTEGRATED SURFACE- AND GROUND- WATER MODEL 4.1
4.2
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.1.1
SWAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.1.2
Integration of SWAT and MODFLOW . . . . . . . . . . . . . . 106
4.1.3
Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2.1
4.2.2
4.3
103
Mathematical description of the soil water module . . . . . . 109 4.2.1.1
Richards equation . . . . . . . . . . . . . . . . . . . . 109
4.2.1.2
Soil constitutive functions . . . . . . . . . . . . . . . . 111
Numerical implementation of SWM . . . . . . . . . . . . . . . 112 4.2.2.1
Soil profile . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.2.2.2
Source and sink . . . . . . . . . . . . . . . . . . . . . . 114
4.2.2.3
Top boundary condition . . . . . . . . . . . . . . . . . 115
4.2.2.4
Bottom boundary condition . . . . . . . . . . . . . . 120
4.2.2.5
Time discretization in SWM . . . . . . . . . . . . . . 121
4.2.3
Surface runoff . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.2.4
Irrigation requirement . . . . . . . . . . . . . . . . . . . . . . . 123
4.2.5
Model interfacing . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.2.5.1
Temporal and spatial discretization . . . . . . . . . . 124
4.2.5.2
Initial conditions . . . . . . . . . . . . . . . . . . . . . 126
4.2.5.3
Input and output files . . . . . . . . . . . . . . . . . . 126
Application of ISGWM . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.3.1
Study Area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.3.2
Model Description . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.4
Simulation results of ISGWM . . . . . . . . . . . . . . . . . . . . . . . 133
4.5
Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 137
x 4.5.1
Innovations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.5.2
Future developments . . . . . . . . . . . . . . . . . . . . . . . . 139
A STREAM DEPLETION RESULTS OF THE LPR MODEL
140
B DATA INPUT INSTRUCTIONS FOR MODFLOW-SDA
153
C DATA INPUT INSTRUCTIONS FOR ISGWM
156
D R SCRIPTS TO GENERATE SWM AND SUW INPUTS
158
Bibliography
162
xi
List of Figures 2.1
Location of study area.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2
Distribution of the land use and land cover in the study area. . . . . . . 12
2.3
Monthly precipitation at meteorological stations. . . . . . . . . . . . . . . 13
2.4
Locations of control points for interpolation. . . . . . . . . . . . . . . . . 20
2.5
Variograms of the elevation of each hydrostratigraphic unit. . . . . . . . 21
2.6
East-west (row direction) cross-sections of the aquifer. . . . . . . . . . . 22
2.7
North-south (column direction) cross-sections of the aquifer. . . . . . . . 23
2.8
Distribution of observation wells. . . . . . . . . . . . . . . . . . . . . . . . 25
2.9
Distribution of registered groundwater wells in the study area. . . . . . 26
2.10 A clip of the map of net corn crop irrigation requirement. . . . . . . . . 28 2.11 Stream network and gauge stations. . . . . . . . . . . . . . . . . . . . . . 31 2.12 River stages along the main course of the Lower Platte River . . . . . . . 32 2.13 Plan view of grid discretization in the Lower Platte River model. . . . . 34 2.14 Time step division in a 30-day stress period. . . . . . . . . . . . . . . . . 35 2.15 Contour map of initial heads and the variogram model. . . . . . . . . . 37 2.16 Distribution of hydraulic property zones in Layer 1. . . . . . . . . . . . . 38 2.17 Distribution of hydraulic property zones in Layer 2. . . . . . . . . . . . . 38 2.18 Distribution of hydraulic property zones in Layer 3. . . . . . . . . . . . . 39
xii 2.19 Distribution of hydraulic property zones in Layer 4. . . . . . . . . . . . . 39 2.20 Distribution of GHB, CHD and RIV grid cells. . . . . . . . . . . . . . . . 42 2.21 Development of groundwater pumping wells between 1950 and 2004. . 43 2.22 Distribution of precipitation recharge zone. . . . . . . . . . . . . . . . . . 44 2.23 Occurrence of groundwater evapotranspiration. . . . . . . . . . . . . . . 45 2.24 MODFLOW packages used in the model. . . . . . . . . . . . . . . . . . . 47 2.25 Average monthly flow rate of each term. . . . . . . . . . . . . . . . . . . 52 2.26 MAE of each observation well. . . . . . . . . . . . . . . . . . . . . . . . . 53 2.27 Flow chart of MODFLOW-SDA. . . . . . . . . . . . . . . . . . . . . . . . 60 2.28 SDF results of the hypothetical model. . . . . . . . . . . . . . . . . . . . . 62 2.29 Distribution of mean errors of SDF results simulated with SDA. . . . . . 67 3.1
Schematic diagram of an example MODFLOW grid superimposed by a stream segment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.2
Plan view showing an example MODFLOW grid superimposed by three connected stream segments. . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.3
Cross-section subdivision for stream depth calculation. . . . . . . . . . . 81
3.4
Algorithm of computation of the intersection between the submerged streambed area and the MODFLOW grid (plan view). . . . . . . . . . . . 83
3.5
Schematic of the hypothetical stream-aquifer model. CS [101] to CS [109] represent the main channel; CS [201] to CS [202] represent the diversion channel; and CS [301] to CS [305] represent the tributary channel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
xiii 3.6
Distribution of stream-aquifer seepage and contours of groundwater hydraulic heads near the streams. Hgw and Qgw denote the groundwater hydraulic head and stream-aquifer seepage, respectively. Negative Qgw represents discharge from the stream to the aquifer and vice versa. . . . 94
3.7
Streamflow and stream stage at each cross-section along the main channel (steady-state simulation with groundwater pumping). . . . . . 96
3.8
Hydrographs simulated with the CSR package at different cross-section in the transient-state model. . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.9
Comparison of simulated outflow hydrographs at CS [109] using CSR, SFR2 and HEC-RAS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.1
Exchange flows between the models. . . . . . . . . . . . . . . . . . . . . . 109
4.2
Uniform movement of overland flow. . . . . . . . . . . . . . . . . . . . . 118
4.3
Top boundary condition types used for SWM. . . . . . . . . . . . . . . . 121
4.4
Interface fractions generated by overlaying MODFLOW grid and SWAT HRUs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
4.5
Location and topography of the Johnson Creek watershed. . . . . . . . . 128
4.6
Land use of the Johnson Creek watershed. . . . . . . . . . . . . . . . . . 129
4.7
Soil distribution and textures of JCW. . . . . . . . . . . . . . . . . . . . . 129
4.8
Calculated and measured θ of soil in JCW. . . . . . . . . . . . . . . . . . 131
4.9
Discretization of SWAT HRUs and MODFLOW grid cells for JCW. . . . 133
4.10 Simulated versus observed monthly streamflows. . . . . . . . . . . . . . 134 4.11 Annual net recharge versus precipitation. The ratio of the annual recharge and precipitation is labeled along the annual recharge line. . . 135 4.12 Mean monthly net recharge distribution. . . . . . . . . . . . . . . . . . . 136 4.13 Mean monthly irrigation pumping volume for different months. . . . . 137
xiv
List of Tables 2.1
Description of hydrostratigraphic units in the modeled area. . . . . . . . 18
2.2
Values of hydraulic property parameters. . . . . . . . . . . . . . . . . . . 49
2.3
Values of precipitation recharge ratios. . . . . . . . . . . . . . . . . . . . . 50
2.4
Cumulative water budget of the LPR model. . . . . . . . . . . . . . . . . 51
2.5
Locations and depletion rates of hypothetical wells . . . . . . . . . . . . 63
3.1
Comparison of MODFLOW stream packages. . . . . . . . . . . . . . . . 76
3.2
Stream channel properties in the hypothetical model. . . . . . . . . . . . 92
3.3
Volumetric water budget of steady-state simulations. Positive and negative values refer to gains and losses to the aquifer, respectively. . . 95
4.1
Values of κ used in overland flow estimation . . . . . . . . . . . . . . . . 119
4.2
Overland runoff qs and time ts to entirely convert to runoff for different ponding depth h0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.3
Major land use types in JCW . . . . . . . . . . . . . . . . . . . . . . . . . 128
4.4
Properties of the first 10 wells with largest pumping rates . . . . . . . . 137
1
Chapter 1 INTRODUCTION 1.1 Background Groundwater had been thought as an isolated resource and was managed separately from the surface water system. Much of the groundwater research performed earlier primarily focused on groundwater. However, groundwater systems are often closely connected to and significantly affected by its upper boundary, the land surface (Sophocleous, 2002; Winter, 1995, 1999). Their interaction constitutes an essential part of the hydrologic cycle. Changes of the natural environment, including the atmosphere and land surface, are affecting and will continue to affect the groundwater systems, directly or indirectly. These impacts on groundwater have become a serious concern for many scientists and managers of water resources. Understanding such interactions between the land surface and groundwater is, therefore, essential for the groundwater management and for ensuring the quality and sustainability of water resources. Traditional methods focusing on either surface water or groundwater are obviously inadequate to address these concerns. In addition, due to their complex nature, it’s difficult to directly observe the in-
2 teraction processes. To overcome these limitations, an Integrated Surface-Ground Water Model (ISGWM), is needed to investigate such processes of the interaction and quantitatively estimate the impacts of their interactions on the surface water and groundwater systems. In most hydrologic models, groundwater components are ignored or oversimplified as a lumped conceptualization described by several individual parameters. Due to such simplifications, compared to surface water resources, sensitivity studies of groundwater systems to the climate and human factors is relatively scarce. These factors include whether conditions, surface water bodies, land cover and land use (LCLU), groundwater pumping, crop planting and irrigation practices. The main forces driving the hydrologic cycle can also have indirect impacts on groundwater systems. Precipitation is the primary component of groundwater recharge. Temperature and radiation are important factors to determine evapotranspiration. The stream stage is an indicator of the interaction between stream and ground water. Change of land cover and land use alters the physical processes on land surface. Vegetation has accordingly different leaf systems and root systems which will impact water consumption demands. Groundwater pumping could directly decrease the groundwater storage significantly. Analysis of the impact of these factors on the groundwater system remains unclear but important for sustainable conjunctive water management. Before the water reaches the groundwater table, it undergoes a number of processes in the hydrologic cycle, such as canopy interception, runoff generation, and movement in unsaturated zones, etc. These processes are essentially complex, therefore making it difficult to account for the effect of all these processes on groundwater recharge. Groundwater recharge increases aquifer storage and raise the groundwater level. In turn, the groundwater level is the bottom boundary for
3 soil water movement. Besides precipitation, irrigation is also a very important groundwater recharge source, particularly in agricultural areas. While part of the irrigation water is used by crops for transpiration or evaporation, part of it infiltrate into the subsurface and eventually become groundwater recharge. However, the pumped volume was not measured in most areas. Thus, estimation of the groundwater pumpage is of great importance for the groundwater recharge model. Uncertainty of groundwater recharge estimation is still high without proper consideration of changes in the groundwater levels and return flow of irrigation water. Water scarcity is a critical issue in arid and semi-arid areas, especially under the threat of increasing water demands and global change (Wheater et al., 2007). In these areas, groundwater, as the primary water supply resource, is being pumped unsustainably (Anderson and Lloyd H. Woosley, 2005). Indicated by McGuire et al. (2003), water levels had declined more than 100 feet in some parts of the High Plains aquifer by 1980 and continued to decline for additional 40 feet in a large area by 1999. In the 2000 water use report (Hutson et al., 2004), Nebraska was estimated to have the largest ground-water withdrawals and 87% of the population depends upon groundwater for water needs, as its semi-arid and arid climate reduces the availability of surface water. According to the Farm and Ranch Irrigation Survey 2007, Nebraska also has the largest irrigated area, and over 90% of groundwater is used for irrigation. Therefore, in Nebraska, groundwater is the most important natural resource and needs to be managed strategically and sustainably. The study area is located in the Lower Plate River Basin which is a highly productive farmlands that include both irrigated and dryland acreages. The Platte River, serving as a valuable water source for Nebraska, flows from east to west across the state. The Platte River basin is of great environmental and ecological
4 significance. The riparian zone along the Platte River provides habitats to some endangered species such as least tern, piping plover, and whooping crane, as well as home to many other plants and animals. In addition, it’s an important irrigated area in Nebraska. However, the stream flow and groundwater could be potentially depleted due to increasing needs of water. To understand the water system and better manage the limited water resources, an integrated water management model tool is needed. Defined in the Nebraska Ground Water Management and Protection Act, Section 46-656.28, hydrologic connection between groundwater and streams exists where pumping will cause streamflow depletion within 50 years greater than 10% of the pumping rate. In the hydrologically connected area, drilling new wells or putting new acres into irrigation may be restricted. The coverage of the hydrologically connected area can be affected by geologic condition and stream properties. Under certain assumptions, the depletion can be estimated analytically (Hunt, 1999). In the Lower Platte basin, however, the geologic setting and stream conditions are far more complex to apply the analytical solutions. To analyze the groundwater flow and its connection with the Lower Platte River, a numerical model is needed.
1.2 Research summary and objectives In this research, an integrated surface-ground water model (ISGWM) is developed and applied to address problems related to interaction between the land surface and groundwater in the Lower Platte River basin. Along with this goal, two other modeling tools have been developed. To be specific, in this dissertation research, three objectives are tried to achieve insights of the surface water and groundwater
5 interaction questions in this research area. These objectives are described below. Stream depletion analysis in the Lower Platte River Basin The Platte River, serving as a valuable water source for Nebraska, and the Platte River basin is of great environmental and ecological significance. Maintaining long-term sustainability of in-stream flow is one of the crucial water management objectives in Nebraska. To evaluate irrigation effect on stream depletion for the Lower Platte river (LPR) basin, a regional three-dimensional, transient groundwater-flow model is established with MODFLOW. Conventionally, calculation of stream depletions requires the performance of a change analysis, both a baseline run and a scenario run of MODFLOW, to calculate the baseflow difference. The process for obtaining a stream depletion distribution pattern is computationally inefficient because separate model runs with a pumping well in each individual grid cell are required. This study proposes a simplified and more efficient solution to estimate a stream depletion rate at a given location by re-using coefficients of the baseline run in the scenario runs. The stream depletion analysis (SDA) tool, therefore linearizes the flow equation because the head coefficients become constant. Since flow coefficients have been calculated in the baseline model run, the solution becomes very efficient and does not require the MODFLOW input files to be repeatedly read into the scenario run. The tool has been successfully applied to the LPR model. A modeling tool to simulate wide-channel streams in MODFLOW The second modeling tool is a cross-section based streamflow routing (CSR) package for MODFLOW to simulate the streamflow and the interaction between streams and aquifers for the stream with a width larger than the
6 MODFLOW grid size. In the CSR package, streams are divided into stream segments which are formed by two consecutive cross-sections. A crosssection is described by a number of streambed points that determine the geometry and hydraulic properties of the streambed. The stream water depth and streamflow at the cross-sections are related by the Single Channel method, the Divided Channel method, a data table or a power function. A rapid algorithm is used to compute the submerged area of the MODFLOW grid. The streambed conductance of a grid cell is computed based on its submerged area, streambed hydraulic conductivity and thickness. Streamaquifer seepage is subsequently estimated as the product of the streambed conductance and difference between the stream stage and groundwater hydraulic head. Stream-aquifer seepage is treated as lateral flow in the streamflow routing computation with the Muskingum-Cunge method or mass conservation method. Develop an integrated model Traditional land surface models oversimplify subsurface water processes and ignore the lateral groundwater flow. Groundwater flow models, on the other hand, cannot depict the physical mechanisms of land surface processes which are often simplified as a constant parameter and need to be externally specified. The third part of this dissertation aims to quantify the impacts of natural processes and human activities on ground-water dynamics in a highly agricultural area by developing an integrated surface-ground water model (ISGWM). ISGWM is developed by coupling SWAT and MODFLOW, which are linked by a soil water module (SWM), developed based on a non-iterative solution to the 1D Richards equation. SWM explicitly represents infiltration,
7 soil evaporation, unsaturated water flow, root water update, and lateral drainage. Water-exchanges across the surface-subsurface and unsaturatedsaturated zone interfaces are defined as the system’s dependent top and bottom boundaries of the soil profile, respectively. Stream-aquifer interaction is simulated by the SFR package in MODFLOW. Taking advantage of the simulation capacities of these models, ISGWM can simulate the physical hydrologic processes in three domains and their interactions. The model can also quantify different components in the hydrologic cycle, such as evaporation, transpiration, and groundwater recharge. Irrigation requirement is also simulated using the soil water deficit or crop water stress.
1.3 Agenda of the dissertation The dissertation contains four chapters. Chapter 1 is the overall introduction providing the background and objectives and summary of this dissertation research. The second chapter presents the first research project which is stream depletion analysis in the Lower Platte River basin. Chapter 3 introduces the cross-section based streamflow routing package for MODFLOW. Chapter 4 demonstrates the development of ISGWM and its application in a small watershed of the Lower Platte River basin. Chapter 3 has been published as a research article in Environmental Modeling and Software (Ou et al., 2013).
8
Chapter 2 STREAM DEPLETION ANALYSIS IN THE LOWER PLATTE RIVER BASIN 2.1 Introduction Groundwater pumping alters the natural water flow dynamics between the surface water and groundwater systems by intercepting groundwater flow originally discharging into surface water body, or by increasing the rate of water movement from the surface water body into an aquifer. Streamflow depletion has become an important concern in Nebraska because extensive use of groundwater may threaten river habitats. An investigation to the effect of groundwater pumping on streamflow is needed. In Nebraska, the Platte River is of great environmental and ecological importance. The Lower Platte River is defined as the section of the Platte River from its confluence with the Loup River to its confluence with the Missouri River. The Lower Platte River basin located in the east-central Nebraska
9 is an area of intensive groundwater development for agricultural irrigation uses. In most part of this area, the surface water and groundwater systems are probably hydrologically connected. Thus, it is necessary to evaluate the potential effects of groundwater extraction on streamflow for water resources management and prediction. A regional, three-dimensional, transient groundwater-flow model can be used as a tool to analyze stream-aquifer interactions induced by groundwater pumping.
2.1.1
Purpose
The purpose of this research is to describe the development and application of the groundwater flow model and a stream depletion analysis (SDA) tool, to analyze the groundwater flow and its connection with the Lower Platte River. The calibrated model is used to evaluate different potential effects of groundwater withdrawals on streamflow in the Platte using these estimated parameters under different hypothetical groundwater extraction scenarios. And SDA is used to defined hydrologically connectedness between the stream and aquifer.
2.1.2
Study area
The study area (model domain) is located in the Lower Platte River basin of eastcentral Nebraska and encompasses portions of Platte, Colfax, Dodge, Washington, Polk, Bulter, Saunders, Douglas and Cass counties, covering an area of 1,930 mi2 (Figure 2.1). The northern and southern boundaries of the model were selected partly to coincide with topographic divides that separate drainage basins. The biggest cities in the area are Columbus and Fremont with the populations of 20,971 and 25,174, respectively, according to the records of U.S. Census 2000. The
10 Platte River flows mainly eastwards between Columbus and Fremont and then southeastwards after Fremont. The lower Platte River ecosystem provides riverine habitat for forage, reproduction, or living space of resident and migratory fish and wildlife species including state- and federally listed endangered species, that is, pallid sturgeon, piping plover, interior least tern, and river otter (Ginting et al., 2008). Other important waterways include the Shell Creek, Salt Creek and Wahoo Creek. The surface elevation decreases west to east from 1718 to 1039 feet above the datum of NAD 1983. Groundwater flow within the model area generally follows the flow of the Platte River, with localized flow occurring towards rivers.
Figure 2.1: Location of study area.
11 Based on the 2005 Land Use Map (Figure 2.2) compiled by the Center for Advanced Land Management Information Technologies (CALMIT) at University of Nebraska-Lincoln, the primary land uses in the Lower Platte River basin consists of dry cropland, irrigated cropland, pastureland, rangeland, and riparian forest and woodlands. The corn and soybean are the major crop types. The dry cropland, pastureland and rangeland are distributed on the highlands, where the groundwater table is usually deeper, while the irrigated cropland is distributed mainly in the Platte River alluvium, where groundwater depth is relatively shallow. Along the main course of the Platte River distribute the riparian forest and woodlands that can intercept the groundwater discharge into the river. There is no large metropolitan area, but significant well fields exist near the Platte River in the study area for Fremont, Omaha and Lincoln. The area possesses a humid continental climate with cold but relatively dry winters and hot and occasionally humid summers. Concurrence of heat and rain in the summer provides ideal growing conditions for crops in this area. The area receives an average annual precipitation of 28.8 inches according to the records at five meteorological measure stations in the cities from west to east of Columbus, Clarkson, David City, Fremont and Wahoo (Figure 2.3). Major precipitation, which is nearly 65%, occurs between May and September. A slight decrease in the annual precipitation can be observed from southeast to northwest, but the average precipitation at each station is similar. In the study area, the principal aquifer consists of saturated unconsolidated sediments of the Quaternary deposits and alluvial aquifers, and the Tertiary Ogallala Group (Nebraska Department of Natural Resources, 2005). The Quaternary deposits cover almost the entire study area and consist of sand, gravel, silt, and clay. Most of the Quaternary deposits are unconsolidated due to the occurrence of
12
Figure 2.2: Distribution of the land use and land cover in the study area. erosion and soil formation during the Quaternary time. Most of the groundwater irrigation wells have screens across this group and large amount of water are withdrawn during the irrigation season. The Ogallala Group mostly consists of sand and gravel with a few silt and clay, and limestone and sandstone sediments at some locations. The Ogallala Group underlies the Quaternary deposits, and the sediments are partially consolidated due to compaction and cementation. Moreover, the alluvial deposits occur in the Platte River valley and some other river valleys in the study area. Similar to the Ogallala Group, the alluvial deposits consist mostly of sand and gravel with a few silt and clay sediments. The alluvial deposits were, however, mostly deposited during the Quaternary time so that they
13 5 C O L U M C L A R K D A V ID F R E M O W A H O
M o n th ly P r e c ip ta tio n ( in )
4 3
O
B U S S O N C IT Y N T
2 1 0
J a n
F e b
M a r
A p r
M a y
J u n
J u l
M o n th
A u g
S e p
O c t
N o v
D e c
A v e ra g e
Figure 2.3: Monthly precipitation at meteorological stations. are unconsolidated. In contrast to the Ogallala Group, the hydraulic conductivity of the alluvial deposits is generally higher due to less compaction and coarser textures. The saturated thickness of the principal aquifer varies from 0 to over 500 feet. In the study area, undifferentiated shale is considered the bedrock unit, which contains thicker shale and thinner shaley chalk. The unit directly underlies the Ogallala Group and the Quaternary deposits. Because the hydraulic conductivity of this unit is generally very small, it is not in hydrologic connection with the Platte River and other rivers in the study area and it is served as a confining layer.
2.1.3
Previous studies
There is a vast volume of studies on streamflow depletion. Wen and Chen (2006) applied nonparametric techniques to analyze the impacts of groundwater irrigation on the streamflow of the major rivers in Nebraska. They found that the decreasing trends of streamflow in these rivers match the pattern of increasing irrigation wells and decreasing groundwater levels. Furthermore, they noted that streamflow depletion in the Republican River Basin may be attributed to that a
14 large amount of groundwater was withdrawn for irrigation uses in the adjacent states of Nebraska, Kansas, and Colorado. Zume and Tarhule (2007) used Visual MODFLOW to investigate the effects of groundwater withdrawal on streamflow of the Beaver-North Canadian River in semiarid northwestern Oklahoma. They divided the basin-wide stream into several segments and noted that the whole study area could be regarded as a gaining system except for one segment. On the other hand, stream-aquifer interactions around an island where water for irrigation use was extracted from streams have been studied by Rodr´ıguez et al. (2006). Their results suggested that a reduction in irrigation water allocations and a switch to more efficient irrigation practices was needed to reduce excessive water table rises on the island. A variety of methods have been reported for the analysis of streamflow depletion, which include analytical solutions developed by Theis (1941), Glover and Balmer (1954), Hantush (1967), Hunt (1999), Butler Jr. et al. (2001), Lough and Hunt (2006), and Sun and Zhan (2007)). These analytical solutions can provide a general idea of the effects of groundwater pumping on streamflow. In current research, numerical models have been used to quantify streamflow depletion instead of analytical solutions (Sophocleous et al., 1995; Chen and Shu, 2002; Nyholm et al., 2002; Chen and Chen, 2003; Hunt and Scott, 2005; Rodr´ıguez et al., 2006; Zume and Tarhule, 2007; Sanz et al., 2011). Numerical models can provide more flexibility in dealing with real hydrologic conditions of actual stream-aquifer-well-riparian vegetation-precipitation systems, such as groundwater evapotranspiration (ET), aquifer heterogeneity and anisotropy, partial penetration streams, heterogeneous hydraulic properties of streambed sediments, and existence of multiple streams. MODFLOW, a widely used groundwater flow model, is usually employed to evaluate the impacts of groundwater exploitation on baseflow variations. River package and Stream package embedded in Visual MODFLOW are two alternative
15 packages to quantify streamflow depletion (Sophocleous et al., 1995; Rodr´ıguez et al., 2006; Zume and Tarhule, 2007). Irrigation is the largest consumer of groundwater use, and most of the groundwater irrigation wells occur in the Platte River valley area. Peckenpaugh and Dugan (1983) developed a two-dimensional groundwater flow model to simulate the impacts of groundwater irrigation in the Central Platte and Lower Loup Natural Resources districts of Nebraska. They examined three management alternatives and found that water-level would decline even without additional groundwater development. However, they did not include the streambed characteristics in their groundwater flow model, and they also assumed that the entire groundwater system was an isotropic unconfined aquifer. Peterson et al. (2008) developed a regional groundwater flow model to simulate the effects of groundwater irrigation on baseflow in the Elkhorn and Loup River Basins in Nebraska. They noted that these effects were minimal before 1970 but increased steadily after 1970. It may be a result of the fact that groundwater irrigation became more common throughout the entire study area after 1970 than the limited areas near the Platte River before 1970. Aquifer horizontal and vertical hydraulic conductivity (aquifer anisotropy) as well as aquifer heterogeneity are of importance to control stream-aquifer interactions (Sophocleous et al., 1995; Chen and Yin, 1999)). Furthermore, streambed hydraulic conductivity has been shown to be the most important hydraulic parameter that affects the stream-aquifer interactions (Sophocleous et al., 1995; Chen and Shu, 2002). Chen et al. (2008) used direct-push technique, conducted by Geoprobe, to collect the electrical conductivity (EC) logs and cores of channel sediments in the Platte River. The EC logs and sediment cores indicated that the streambed sediments below the channel surface consist mostly of sand and gravel, thus the Platte River was considered hydrologically connected with adjacent aquifers. Wen and Chen
16 (2006) pointed out that 4 gauge stations have decreasing trends of streamflow over a total of 27 gauging stations in the Platte River during the period of 1950 to 2003, while other stations have insignificant trends. Consequently, a more accurate streamflow analysis in the Lower Platte River basin is beneficial to understand the stream-aquifer interactions with new findings derived from a groundwater flow model.
2.2 Conceptual groundwater-flow model This section provides information on the conceptual groundwater flow model. The conceptual model of the Lower Platte River basin, upon which the numerical model is based, is a simplified description of the aquifer system. Data were compiled to develop the conceptual groundwater model. Data were retrieved from the U.S. Geologic Survey (USGS), Nebraska Department of Natural Resources (DNR) and other state agencies and organizations.
2.2.1
Hydrostratigraphic units
2.2.1.1
Hydrogeologic settings
The geologic setting of the area is described in the Study area section. Based on the geologic setting, the studied aquifer system is divided into four hydrostratigraphic units (Table 2.1). The first layer includes the top soil layer and is composed mainly of silt and clay sediments, and alluvial sand and gravel deposits in the Platte River valley and Quaternary sand and gravel sediments in the upland area, with a few till sediments embedded. The second layer consists of silt, clay and fine sand sediments, which is regarded as a confining layer, while the third layer is
17 mostly composed of Ogallala sand and gravel sediments as well as the Dakota sandstones. The bottom layer mostly contains silt and clay with a few limestone and sandstone sediments. The sediments below the bottom layer are considered as bedrock unit and are not hydrologically connected with the Platte River; therefore, these sediments were not included in the model.
Table 2.1: Description of hydrostratigraphic units in the modeled area.a System
Hydrologic unit
Composition descriptions
Quaternary
Top soil Platte River alluvium Loess
silt and clay with sandy deposits alluvial sand, gravel and silt deposited within incised bedrock valley of the Platte River silt with a small amount of very fine sand and clay deposited as wind-blown dust glacial deposited silty, sandy clay with some gravel, pebble, and cobbles
Till Tertiary Cretaceous
Dakota group
gravel, sand, silt, clay, with some limecemented beds massive to crossbedded friable sandstone interbedded with clayey to slightly sandy shales; sandstone may contain ironstone or siderite concretions and chert pebbles
Maximum thickness (feet)
Layer 1
355
Layer 2
109
Layer 3
148
Layer 4 Undifferentiated shale, shales interbedded with limestone and sandlimestone and stone; shales may be laminated clayey, sandy, sandstone calcareous, fossiliferous and may contain gypsum; limestone may be massive, geodal, fossiliferous or contain chert; sandstone is generally thin bedded and may contain coal. modified from Nebraska Department of Natural Resources (Nebraska Department of Natural Resources, 2005).
Permain and Pennsylvanian Undifferentiated
a
Ogallala group
Hydrostratigraphic unit
109
18
19 2.2.1.2
Elevation calculation
The elevation of the hydrostratigraphic units is determined through the Kriging interpolation method. The procedure is described as follows: 1. Geologs of the registered wells from NDNR and lithology of the state test holes from CSD in the study area were selected as control points (Figure 2.4). The geologs and lithologic records give the depth and the descriptions of the sediment unit at the registered well or test hole location. The elevation of each layer at each control point was determined based on grouping the similar and adjacent sediment composition units. The surface elevation at the USGS groundwater observation wells was also used when calculating the top elevation of the first layer. 2. The spatial variability was analyzed through fitting the semivariogram of the synthesized elevation results with theoretical variogram models. For a two-dimensional dataset, the semivariogram is: 1 γ(h) = 2N (h)
N (h)
∑ [z(u) − z(u + 1)]2
(2.1)
1
where γ(h) is the semivariogram; h is the distance interval along all or some specified directions; N (h) is the number of pairs of data points with such distance interval h; and z(u) and z(u + h) are a pair of values with the distance interval h. 3. The Kriging algorithm was used to interpolate the values at control points according to the fitted theoretical variogram model. Kriging is widely used in hydrogeology as a preferred method to construct gridded datasets (Tonkin and Larson, 2002). It is a minimum variance estimation algorithm, and exact
20
Figure 2.4: Locations of control points for interpolating the elevation of hydrostratigraphic units. interpolator at control points, that is, the interpolated values are equal to the measured values at the locations of the measure points (Kitanidis, 1997). The fitted variograms for the top and bottom elevations of Layer 1, and bottom elevations of Layer 2, 3 and 4 are shown in Figure 2.5. The black dot line is the variogram directly calculated with Equation 2.1 based on the elevation at each control points (or experimental variogram). The blue solid line is the fitted theoretical variogram model. The nugget effect, which may be caused by measurement errors or undetected microvariabilities (Kitanidis, 1997), was found for the bottom models of all layers. The equations of fitted theoretical models are also shown in Figure
21
Figure 2.5: Variograms of the elevation of each hydrostratigraphic unit. 5. The models of (a), (b) and (c) in Figure 5 belong to the nonstationary power model, while (d) and (e), detrended before the variogram calculation, belong to the stationary exponential model. Figure 2.6 and Figure 2.7 display the cross-sections of the hydrostratigraphic units along the row (east-west) and column (north-south) directions. Overall, the elevations of all the layers tend to coincide with the change of land surface, that is, the elevation decreases from west to east. This topographic characteristic determines the overall groundwater flow direction in this area. The Platte River Valley has the lowest elevation, indicating groundwater flows toward and recharge the river. Layer 1 and Layer 3 are relatively thick compared to the other layers. They are the primary groundwater storages.
22
Figure 2.6: East-west (row direction) cross-sections of the aquifer.
23
Figure 2.7: North-south (column direction) cross-sections of the aquifer.
24
2.2.2
Groundwater levels
According to the Nebraska Statewide Groundwater Level Database, there are 572 monitoring wells within the study area (http://snr.unl.edu/data/water/ NebGW_Levels.asp). The wells with less than five observation records are excluded from calibration. For multiple wells located in the same grid cell, only the well having the most observation records were used in the model calibration. Totally 189 wells with 14,153 observation records were selected for the model calibration based on their locations and numbers of observation records. Figure 2.8 displays the distribution of the observation wells and the topographic regions. There are 135 wells located in the Platte River Valley which is a developed agricultural area with densely distributed pumping wells. The change of groundwater levels is caused by the stresses such as recharge, evapotranspiration or groundwater pumping. In this area, the average slope of the linear trends of the groundwater level changes at all the observation wells is 0.098 ft/year, that is, an average increase of 0.4% between 1950 and 2004. This suggests that there was no significant trend of the groundwater level change in this area over the 55 years from 1950 to 2004. The intra-annual variation caused by the irrigation pumpage, however, is significant, and the range of the variation can reach over 60 feet (at Well 128).
25
Figure 2.8: Distribution of observation wells.
26
Figure 2.9: Distribution of registered groundwater wells in the study area.
2.2.3
Groundwater pumping
Groundwater pumping causes the reduction of groundwater resources. In the study area, groundwater pumpage is the largest discharge from the groundwater system. Excluding those wells with a pumpage less than 50 GPM, there are over 5000 registered pumping wells in the study area and over 95% of them are agricultural irrigation wells (Figure 2.9). Therefore, the irrigation pumpage is the primary stress inducing the variation of the groundwater system. However, the pumping volume was not measured in the study area during the modeling period. Thus, estimation of the groundwater pumpage is of great importance for the groundwater model. We estimated the groundwater pumpage based on the irrigation requirement and the pumping capacity of each well.
27 The irrigated area and the water pumping discharge capacity of each well can be obtained from the DNRs registered well database. A ratio of the actual pumpage and the capacity could be calculated if the irrigation requirement is known during the crop growing seasons. The amount of groundwater pumped from the aquifer was assumed to be equal to the irrigation requirement, thus:
β=
αR ∑1n Ai ∑1n (Ci T )
(2.2)
where β is the pumping ratio of the actual pumpage and pumping capacity in the study area; α is a unit conversion factor which is equal to 18.86; R is the annual irrigation requirement in this area, (inch); n is the total number of irrigation wells in this area; Ci is the average pumping capacity of the ith wells in this area (GPM); T is the time of the irrigation period which we assumed to be 92 days in June, July and August every year; and Ai is the irrigated acres of the ith well in this area (acre). Therefore, the pumping rate of a irrigation well is βCi . If considering the total precipitation P (inch) occurring during the irrigated months, we can calculate the pumping ratio as follows: α( R − P) ∑1n Ai β= ∑1n (Ci T )
(2.3)
β was estimated to be 0.183 using Equation 2.3 based on the properties of each well in the study area and the annual net irrigation requirement shown in Figure 2.10, which was compiled by DNR (http://dnr.ne.gov/SurfaceWater/ CountyMapIrrReq.pdf). An average value of 0.13 for β was used after accounting for other crops having less irrigation requirement.
28
Figure 2.10: A clip of the map of net corn crop irrigation requirement in Nebraska (inch/year).
2.2.4
Recharge
Before the water reaches the groundwater table, it undergoes a number of processes in the hydrologic cycle, such as canopy interception, runoff generation, and movement in unsaturated zones, etc. These processes are essentially complex, therefore making it difficult to account for the effect of all these processes on groundwater recharge. A relatively simplified method was used to quantify the groundwater recharge: Q f = λF
(2.4)
where Q f is the flux of groundwater recharge; λ is a ratio specifying the portion of water reaches the water table and becomes recharge to the groundwater system; and F is the rate of water applied to the land surface. F is comprised mainly of
29 precipitation and irrigation, thus F = P + I; P denotes precipitation, and I denotes the amount of irrigation water, that is, the groundwater pumpage for irrigation wells. However, while precipitation occurs on the whole area, irrigation occurs only on the farms during the crop growing seasons. Thus, the two components are estimated separately and subsequently summed up.
2.2.5
Evapotranspiration
In terms of groundwater evapotranspiration (ET), most previous studies neglected it in their simulations since the effects of ET were usually minimal in a regional model (Zume and Tarhule, 2007). Groundwater ET was often specified to occur in the riparian vegetation areas. Chen and Young (2006) estimated that the groundwater ET rate ranged from 0.78 to 1.0 in/yr for the riparian vegetation and that the extinction depth for ET was 15 ft in the central Platte River Valley, Nebraska. Peterson et al. (2008) used a uniform evapotranspiration rate of 14 in/yr and an extinction depth of 5 ft in their model of the Elkhorn and Loup River Basins. In our study, the groundwater ET losses mainly come from the riparian vegetation and occur in the Platte River riparian valley area. We assumed that the potential groundwater ET rate was 27 in/yr in the riparian vegetation zones, and that no groundwater was lost through ET when the depth to water table is greater than 15 feet, that is, the ET extinction depth.
2.2.6
Surface water
2.2.6.1
Stream network
The stream network was generated through delineating a 30-meter USGS DEM in a larger area covering the whole Lower Platte River basin. An ArcGIS extension,
30 ArcSWAT, which originally is a graphical user interface for the SWAT model (Arnold et al., 1998), was used to delineate the watershed and generate the stream network (Figure 2.11). Streams were divided into reaches based on the stream length and accumulation area. The stream network is composed of 34 and 172 reaches of the Platte River and its tributaries, respectively. Each reach has the same channel width and depth which are calculated by ArcSWAT using the following equations: W = 7.49A0.6
(2.5)
D = 0.624A0.4
(2.6)
where W is the channel width of a reach (ft); A is the accumulation area of the reach (mi2 ); and D is the depth of the reach (ft). The calculated widths of the streams were subsequently randomly selected to compare with the visual measurement on the aerial photograph. We found that the calculated results were close to the visual measurements. The streambed conductivity was assigned to each reach based on our field tests (Chen et al., 2008), such that the conductance could be calculated as (McDonald and Harbaugh, 1988): CRIV = Kb LW/M
(2.7)
where CRIV is the conductance of streambed (ft2 /day); Kb is the streambed conductivity (ft/day); L is the stream length (ft); W is the stream channel width (ft); and M is the streambed thickness (ft). We assumed M = 3 ft in this model. The exchange between the stream and aquifer can then be computed on the basis of the conductance and the difference of the stream stage and groundwater level.
31
Figure 2.11: Stream network and gauge stations. 2.2.6.2
Stream stages
Due to the data availability, the river stages along the main course and tributaries of the Platte River are handled by different methods. There are five hydrologic gauge stations along the main course of the Platte River in this area (Figure 2.11). In each stress period, the river stage between gauge stations is linearly interpolated based on the observed water levels at the adjacent upstream and downstream gauge stations. The calculation equation is as follows:
Hr = Hr1 +
L1 ( Hr2 − Hr1 ) L1 + L2
(2.8)
32
Figure 2.12: River stages along the main course of the Lower Platte River in Dec, 2004. where Hr is the river stage at a grid cell; Hr1 and Hr2 are river stages at the adjacent upstream gauge and downstream gauge, respectively; and L1 and L2 are the lengths of river sections from the calculated river grid cell to the adjacent upstream gauge and downstream gauge, respectively. Figure 2.12 demonstrates the river stages along the main course of the Lower Platte River in Dec, 2004. The stages at five gauge stations are nearly linearly distributed with the corresponding downstream river lengths from the inlet point of the Platte River at the boundary of the study area. The overall slope of the river stage is 0.9%. Observed river stages are not available for most tributaries in the study area. Two assumptions were accordingly made: 1. the average annual river stage was consistent with the average annual groundwater table contiguous to the river grid;
33 2. the relative seasonal fluctuation was consistent with that of the nearest gauge station along the main course of the Lower Platte River. With such assumptions, the river stages along the tributaries can be calculated.
2.3 Development of the numerical model The three-dimensional finite-difference modular groundwater-flow model code MODFLOW-2005 (Harbaugh, 2005) was used to simulate the groundwater flow system and the interaction between the Platte River and aquifers. The Lower Platte River groundwater-flow model (LPR) was constructed based on the conceptual model described in the previous section for a period from 1950 through 2004. The data were translated and assigned to the model at discrete intervals in space and time with the assistance of the Visual MODFLOW and Microsoft Excel.
2.3.1
Space and time discretization
MODFLOW employs a finite-difference method to solve the partial-differential equation governing the groundwater flow. The method discretizes the space into a grid framework composed of rectilinear blocks called cells; and divides simulation time into stress periods, in which the external stresses are constant. The length and time units of the LPR model are foot and day, respectively. The LPR model was horizontally subdivided into 180 rows by 201 columns, and vertically into 4 layers. The cell sizes are 1/8 mi by 1/8 mi near the Platte River and gradually increase up to 1/2 mi by 1/2 mi as the distance to the river increases (Figure 2.13). Finer grid was employed near the Platte River based on such reasons: (1) importance of the Platte River in our study; (2) stress from the
34 variation of the river stage; and (3) higher accuracy. The width of the narrowest part of the Platte River on the aerial photograph is approximately 660 ft, which is similar to the smallest horizontal spacing of the grid cells. Mehl and Hill (2010) suggested that refined grids could be used to improve the simulation accuracy when modeling the flux between rivers and aquifers. The interpolated elevations of each model layer were assigned to the center of each grid cell.
Figure 2.13: Plan view of grid discretization in the Lower Platte River model. The model is composed of four layers. To account for the change of confined and unconfined conditions at some location where hydraulic head significantly drops, we used Type 3 aquifer layer (McDonald and Harbaugh, 1988) to simulate each model layer. In MODFLOW, for Type 3 aquifer layer, “transmissivity of the layer varies and is calculated from the saturated thickness and hydraulic conductivity.
35 The storage coefficient may alternate between confined and unconfined values. Vertical flow from above is limited if the aquifer desaturates” (Harbaugh, 2005). At each well location, the status of the grid cell was carefully examined to assure that the gird cell is not “dry” during pumping seasons. The model was simulated for transient conditions and the time was divided into discrete time intervals called stress periods. The simulation period extends from January 1, 1950 to December 31, 2004, for a total 55 years and 660 monthly stress periods. Simulated hydrologic stresses, such as recharge, pumpage and river stages, were updated between stress periods. The hydraulic properties, such as conductivities, specific storages and specific yields, are kept constant during the whole simulation. Each stress period, in turn, was subdivided into a number of calculation time steps. In the LPR model, there are 10 time steps in each stress period. The time steps within each stress period increase with a geometric ratio, for example, 1.2 in our model. Figure 2.14 shows the division of time steps in a 30-day monthly stress period. The length of the early time step is shorter to achieve numerical stability, because the updated external stress at the start of a new stress period can cause rapid change of the hydraulic heads. With the increase of elapsed time, the rate of change in head typically decreases, so time steps can be increased. 1
2
3
4
0 1.2 2.5 4.2
5 6.2
6 8.6
7 11.5
8 14.9
9 19.1
Time Step
10 24.0
30.0
Figure 2.14: Time step division in a 30-day stress period.
Day
36
2.3.2
Initial condition
Initial conditions are important for transient models. In this model, the distribution of initial hydraulic heads was generated through the same interpolation method as described in the previous section to estimate the elevations of hydrostratigraphic units. Because some wells did not have observation records in 1950, the starting year of the model, for the initial hydraulic head at these wells, it needs to be assumed an initial value. The linear regression shows insignificant trend in the groundwater level change, thus for those wells lacking the data in 1950, we used the average measured records between 1950 and 1970; otherwise, for those wells lacking records before 1970, we used the average groundwater level after 1970. Figure 2.15 shows the contour map of the initial head and the fitted variogram model, which is a nonstationary power model.
2.3.3
Hydraulic properties
The required hydraulic properties of aquifers for MODFLOW include horizontal and vertical hydraulic conductivities (Kh and Kv ), specific storages (Ss ) of confined aquifers and specific yields (Sy ) of unconfined aquifers. The horizontal conductivities along the row and column directions were assumed to be identical. Anisotropy occurs only between vertical and horizontal directions. These hydraulic properties define the ability of aquifers to transmit and store groundwater. In reality, these properties change with locations and can be estimated in the field. However, it’s not realistic to measure them at every grid cell. Therefore, these properties were initially estimated based on the values in previous studies, and the sediment properties at the registered pumping wells and the state test holes. Furthermore, two aquifer pumping tests were conducted in the study area and
37
Figure 2.15: Contour map of initial heads and the variogram model. they provided hydraulic property values for the alluvial sand/gravel layers and silt/clay layers. In the model, they were further calibrated to match the hydraulic heads at observation wells. It’s difficult to estimate the hydraulic properties at each individual grid cell. Thus, we divided the model domain into hydraulic property zones (HPZ) in each layer on the basis of the sediment properties and manual adjustment to match the observed hydraulic heads. The hydraulic properties are assumed to be homogeneous within a hydraulic property zone, that is, all grid cells in the zone have the same values of Kh , Kv , Ss and Sy , respectively. During calibration, the values of the properties in each zone were changed. The distribution of the hydraulic zones is shown in Figure 2.16 to Figure 2.19. Through the modeling tests, we found that most wells are more sensitive to
38
Figure 2.16: Distribution of hydraulic property zones in Layer 1.
Figure 2.17: Distribution of hydraulic property zones in Layer 2.
39
Figure 2.18: Distribution of hydraulic property zones in Layer 3.
Figure 2.19: Distribution of hydraulic property zones in Layer 4.
40 the hydraulic properties in Layer 1, partially due to the thickness of the layers. Therefore, to better represent the hydraulic properties in Layer 1, we used four hydraulic property zones for the first layer and two zones respectively for the rest layers as shown in Figure 2.16. Hydraulic Property Zone 1 (HPZ-1) represents the medium to course sand in the Platte River alluvium in Layer 1. HPZ-2 represents the silt, fine sand, sandstone in the Platte River alluvium and highlands in Layer 1. HPZ-3 represents the very high permeable deposits, such as gravel and course sand in the Platte River alluvium in Layer 1. HPZ-4 represents the low permeable deposits, such as till and clay in Layer 1. In the other three lower layers, each layer has two HPZs representing relatively low and high permeable sediment materials, respectively. In Layer 2, HPZ-5 represents the silt, shale and glacial till, while HPZ-6 represents the fine sand and medium-grained sandstone. In Layer 3, HPZ-7 represents the Ogallala sand and gravel while HPZ-8 represents the sandstone or silt of the Dakota group. In Layer 4, HPZ-9 represents the primarily distributed silt, sand and shale deposits, whereas HPZ-10 represents a small portion of limestone and sandstone sediments.
2.3.4
Boundary conditions
Boundary conditions define the locations and flows of water entering and exiting the model domain. Hydraulic head and groundwater flow change as boundary conditions change. There are three types of boundary conditions: (1) specific head such as the constant-head boundary; (2) specific flux such as the recharge and well pumpage; and (3) head-dependent flux such as the evapotranspiration, general-head boundary and river leakage. In terms of the locations, they can also be grouped into: (1) lateral flow along the model boundary; (2) areal vertical
41 recharge and evapotranspiration over the model domain; and (3) internal source and sinks such as rivers and pumping wells. 2.3.4.1
Lateral boundary conditions
Lateral boundary conditions include the constant-head boundary (CHD) and general-head boundary (GHB). CHD specifies the hydraulic head at the grid cells, while GHB calculates the flux at its grid cells by the means of specifying the conductance and the head of an external source (McDonald and Harbaugh, 1988). The hydraulic head can change during the simulation at the general-head grid cells. The flux is the product of the conductance and the difference between the hydraulic head at the grid cell and the head of the external source. We used GHB to represent the whole boundaries except for using CHD for the west part which coincides with the Elkhorn River Figure 2.20. The boundaries of north and south side coincide with the watershed boundaries and are parallel to the overall groundwater-flow direction, while the boundaries of east side are perpendicular to the overall groundwater-flow direction. Thus, we specified low GHB conductance for the northern and southern boundaries, and high GHB conductance for the western boundaries, respectively. In the model domain, our analysis of the groundwater level change during the modeling period suggests an insignificant trend in the water level. Thus, the initial head was used as the specified hydraulic head for GHB and CHD. 2.3.4.2
Groundwater recharge
Similar to hydraulic property zones, we divided the study area into a number of recharge zones to estimate the precipitation contribution to recharge. Each recharge zone was assumed to represent an area with a uniform recharge ratio,
42
Figure 2.20: Distribution of GHB, CHD and RIV grid cells. that is, α in Equation 2.3. However, it is not suitable to use a constant α to represent the infiltration processes throughout the modeling time, because the land use and land cover of each recharge zone had been changing during such a long period, 55 years from 1950 to 2004. Figure 2.21 indicates that the number of the groundwater wells kept growing, implying the increasing development of the agricultural lands in this area (only those wells with a pumpage of 50 GPM or higher are taken into account). Thus, the modeling time was consecutively split into seven phases: 1950∼1955, 1956∼1966, 1967∼1971, 1972∼1977, 1978∼1989, 1989∼1993, and 1994∼2004. As a result, each recharge zone has seven infiltration ratios corresponding to these phases. These infiltration ratios can simplify the infiltration processes in the model, but they are difficult to be estimated directly. Thus, they are one group of the parameters needed to be calibrated manually.
43
Number of registered wells
6000 5000 4000 3000 2000 1000 0 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 Year
Figure 2.21: Development of groundwater pumping wells between 1950 and 2004. According to the Map of land use and land cover in Figure 2.2, the precipitation recharge zones (PRZ) were defined as shown in Figure 2.22. PRZ were initially determined based on the land use and land cover, and adjusted subsequently according to the groundwater levels in observation wells. The Precipitation Recharge Zone 1 (PRZ-1) represents locations of moderate infiltration capacity. It is mainly composed of irrigated farm land with relatively sparsely distributed irrigation wells. PRZ-2 represents the same land cover type as PRZ-1 but with lower infiltration capacity. PRZ-3 and PRZ-4 represents zones with high and very high infiltration capacity, respectively. These two zones represent the riparian zone along the Platte River and the highest developed agricultural areas having densely distributed irrigation wells or irrigation wells of high pumping capacity. PRZ-5 represents the dryland, rangeland and other area of low infiltration capacity. 2.3.4.3
Groundwater evapotranspiration
Evapotranspiration is assumed to occur only at the locations where the groundwater table is shallow. In the LPR model, the groundwater evapotranspiration
44
Figure 2.22: Distribution of precipitation recharge zone. occurs along the Platte River and in the parts of Wahoo Valley (Figure 20). The actual groundwater evapotranspiration was estimated with the evapotranspiration package (EVT) of MODFLOW (McDonald and Harbaugh, 1988). It assumes that evapotranspiration varies linearly with water-table elevation, such that:
ET =
PET,
h − Ze PET , Zs − Ze 0,
h ≥ Zs Ze < h < Zs
(2.9)
h ≤ Ze
where ET is the rate of groundwater ET; PET is the potential groundwater ET; Zs is the ET surface elevation which coincides the land surface in the model; Ze is the extinction elevation below which the groundwater ET ceases; and h is the
45
Figure 2.23: Occurrence of groundwater evapotranspiration. groundwater level. 2.3.4.4
Streambed seepage
Streambed seepage is the product of the conductance and the water level difference between the stream and aquifer. The exchange flux direction is determined based on which of the water levels is higher. If water flows from aquifer to streams, it becomes the baseflow for the stream. Streams also recharge aquifers at the locations where the stream stage is higher than the groundwater level. Decline of the groundwater level induces reduced baseflow or increased stream discharge, causing the streamflow depletion in the streams. The distribution of river cells is shown in Figure 2.20. The Platte River and its tributaries were simulated using the River package of MODFLOW (RIV).
46 2.3.4.5
Pumping wells
Groundwater pumpage was simulated with the Well package of MODFLOW (WEL). The pumping rate was determined based on the annual irrigation requirement and pumping capacity of the well. The pumping period of irrigation wells spans from June to August, that is, three months every year. And the total pumpage was spread into the irrigation seasons. The starting time and ending time of pumpage were determined based on the completion date and abandoned date of the well. And the screen casing depth of the well was used to decide the layer of the pumping occurrence.
2.3.5
MODFLOW Packages
MODFLOW is composed of different packages, which are divided based on their hydrologic representation and computation functions. The packages we used in the model are shown in Figure 2.24. The packages can be grouped into four major types. The first type is the Stress packages. The packages of this type represent the boundary conditions causing the variation of the groundwater system. It includes the General-Head Boundary (GHB) Package, Well (WEL) Package, Evapotranspiration (EVT) Package, Recharge (RCH) Package, River (RIV) Package and Constant-Head (CHD) package. GHB was used to simulate the sources and sinks at the boundary of the model domain. WEL was used to represent the groundwater pumpage for all kinds of uses. EVT and RCH were used to estimate the rate and location of groundwater evapotranspiration and recharge, respectively. RIV was used to simulate the variation of the rivers in the study area, whereas CHD represent the Elkhorn River in the northeast corner.
47
Stress
Internal
Control
OB
GHB WEL
LIST
BAS
EVT
DIS
OC
HYD
RCH
BCF
SIP
HOB
RIV
ZONE
RVOB
CHD
Figure 2.24: MODFLOW packages used in the model. The second type, the Internal packages, describes the properties of the groundwater system and simulates the internal groundwater flow. It includes the Basic (BAS) Package, Discretization (DIS) Package, Block-Centered Flow (BCF) Package and Zone (ZONE) package. BAS reads the initial head and the location of active model domain. DIS is for the discretization file defining the length and time units, and spatial and temporal discretization of the model. ZONE is simply the zone distribution of hydraulic properties and evapotranspiration. The third type, the Control packages, controls the model output format and solves the groundwater flow equation. It includes the Output-Control (OC) Package and Strongly Implicit Procedure (SIP) Package. OC is used to define the output format and location of simulated hydraulic head and overall volumetric water budgets. SIP is a method for solving the groundwater flow equation. The forth type, the Observation packages, outputs the model results according to the specifications. It includes the List (LIST) package, Hydrograph (HYD) package, Hydraulic Head Observation (HOB) package and River Observation
48 (RVOB) package. LIST file contains all the general model information. HYD outputs a hydrograph of the head at specified grid cells. HOB output the calculated head vs. observed head. And RVOB is used to estimate the exchange flux between the river and aquifer. An external program called Zone Budget was also applied to compute the flows in specified zones (Harbaugh, 1990).
2.4 Model calibration The processes of calibration were to change the parameter values to match the model results with observations. The model was calibrated manually through the “trial-and-error” approach. During the manual calibration, the parameter zones and initial parameter values are defined and estimated. Observations used in model calibration include the observed hydraulic heads in the 189 observation wells. Quantitatively, the calibration objective can be expressed as: n
SSQ =
∑ ( Hi − Oi )2
(2.10)
i =1
where SSQ is the sum of square of modeled residuals; Hi is the ith modeled head; and Oi is the ith observed head. The goal of calibration is to minimize the calibration objective, that is, SSQ. Another model evaluation criterion is the mean absolute error (MAE), which can be calculated as: n
MAE =
∑ | Hi − Oi |
i =1
(2.11)
49
2.4.1
Parameters in model calibration
Parameters needed to be calibrated include: (1) parameters of hydraulic property, such as Kh , Kv , Ss and Sy ; (2) precipitation recharge ratios, α. Theoretically, each grid cell possesses its own values for these parameters. It is not, however, feasible either economically or technically to examine the parameter value at each single grid cell. Parameter zoning is commonly used to group the grid cells with the same or similar properties, therefore reducing the difficulty in parameter calibration but still correctly representing the reality to the utmost. The manual trial-and-error process was used for the calibration to identify values of parameters that provided a reasonable fit between simulated and measured hydraulic heads within the model domain. More importantly, the manual calibration defined the parameter zones. The zone division of each parameter is described in Chapter 2.3. Initial parameter values calibrated through the trial-and-error method are shown in Table 2.2 and Table 2.3. The SSQ is 2.2 × 105 and the MAE is 2.76 ft. Table 2.2: Values of hydraulic property parameters estimated through calibration. Zone Number
Kh (ft/day)
Kv (ft/day)
Ss (ft−1 )
Sy
HPZ 1, Layer 1 HPZ 2, Layer 1 HPZ 3, Layer 1 HPZ 4, Layer 1 HPZ 5, Layer 2 HPZ 6, Layer 2 HPZ 7, Layer 3 HPZ 8, Layer 3 HPZ 9, Layer 4 HPZ 10, Layer 4
164 0.3 321 0.02 0.164 16.4 39.37 3.28 0.016 32.8
11.5 0.03 33 0.0016 0.0164 1.31 3.28 0.164 0.001 3.28
3.0 × 10−5 3.0 × 10−4 2.0 × 10−5 1.0 × 10−5 1.0 × 10−5 2.0 × 10−5 2.0 × 10−5 1.0 × 10−5 3.0 × 10−6 3.0 × 10−5
0.22 0.08 0.18 0.06 0.06 0.15 0.20 0.08 0.02 0.10
50 Table 2.3: Values of precipitation recharge ratios estimated through calibration. Zone Number PRZ 1 PRZ 2 PRZ 3 PRZ 4 PRZ 5
50∼55
56∼66
0.05 0.02 0.10 0.18 0.01
0.07 0.02 0.15 0.24 0.01
Year 67∼71 72∼77 0.07 0.05 0.12 0.25 0.01
0.07 0.05 0.21 0.26 0.01
78∼89
89∼93
94∼04
0.09 0.07 0.21 0.26 0.03
0.11 0.07 0.22 0.26 0.03
0.13 0.10 0.22 0.28 0.03
2.5 Model results The calibrated LPR model was evaluated to be accurate compared with the observed groundwater levels. The regional water budget and model fit are discussed in this section. The stream-aquifer connectedness can be evaluated through hypothetical simulations with the calibrated model and the SDA tool.
2.5.1
Regional water budget
Table 3.3 lists the cumulative water budget of the LPR model between 1950 and 2004. In the table, the STORAGE term refers to the water storage change of the aquifer. The percentages in the parentheses are the ratio of the term and total inflow or outflow. We can observe that the positive storage change (inflow) is nearly equal to the negative storage change (outflow), meaning that the total storage of the aquifer did not change significantly. This can also be proved by the relatively stable groundwater levels in the observation wells. In the inflow terms, recharge is the primary water resource for the groundwater system; while the groundwater pumping and river leakage are most important in the discharge terms. The percent discrepancy between the total inflow and outflow is about 0.01%, indicating the numerical solution preserved good precision.
51 Table 2.4: Cumulative water budget of the LPR model. Flow In
Flow Out %
Term
ft3
%
ft3
STORAGE CONSTANT HEAD WELLS RIVER LEAKAGE ET HEAD DEP BOUNDS RECHARGE TOTAL
3.811E+11 7.671E+09 0.000E+00 2.838E+11 0.000E+00 6.587E+10 5.974E+11 1.336E+12
29 1 0 21 0 5 45 100
3.886E+11 3.363E+10 3.980E+11 4.464E+11 4.958E+10 1.958E+10 0.000E+00 1.336E+12
29 3 30 33 4 1 0 100
Figure 2.25 shows the average monthly water flow rate of each term. The terms in the parentheses are outflow from the aquifers and have negative values in the figure. The seasonal variability of recharge and well pumping is outstanding due to the irrigation schedule and the seasonal distribution of rainfall in this area. Inflow and outflow of river leakage are relatively stable due to the small variation of the groundwater level and river stage.
2.5.2
Groundwater levels
The simulated time series of hydraulic head at observation wells are output with the package HYD, while the simulated values for each observation are interpolated at the specific measured time and subsequently output with the package HOB (Figure 2.24). Overall, the mean absolute error (MAE) of the calibrated model for the observed groundwater levels is 2.76 ft. The fit of simulated to observed hydraulic heads is generally good, and we consider the model achieved a reasonable representation of the groundwater flow system in this area. Figure 2.26 shows the spatial distribution and histogram of the MAE of each single well. The area of high uncertainty is delineated. The fit is relatively
52
Figure 2.25: Average monthly flow rate of each term. poor in these parts. The high uncertainty is considered to be attributed to: (1) observation data shortage; (2) geologic data shortage; (3) irrigation data shortage; (4) boundary condition simplification; and (5) difficulty of simulating the wells with variations spanning over a wide range. The uncertainty areas are located close to the model boundary. Thus, the simplification of the boundary conditions, using mean hydraulic heads, will affect the simulation performance, especially for the wells having very rapid temporal variations of the groundwater level to 50 ft such as Well 111. There are some irrigation wells in these areas, but the properties (such as pumping capacity) of the pumping wells are not documented. However, the mean groundwater levels were well simulated at these locations. The MAEs of over 60% and 80% of the observation wells are within 3 ft and 5 ft, respectively. We can observe that the MAE is smaller as the distance from the Platte River is getting shorter. This provides higher confidence in the streamflow depletion analysis for the Platte River.
53
Figure 2.26: MAE of each observation well.
54
2.6 Streamflow depletion One of the most important objectives of this study is to analyze the hydrologic connection between the Platte River and the aquifer through accurate numerical modeling. In Nebraska, hydrologically connected area is defined as the area in which pumping of a well for 50 years will deplete streamflow by at least 10 percent of the pumped amount in accordance with NDNR Rules of Practice and Procedure, Title 457 Neb. Admin. Code Chapter 24 §001.02. Thus, stream depletion factor, which is defined in Equation 2.12, is an indicator for hydrologic connection between the two water systems.
SDF = 100
Qr0 − Qr Qp
(2.12)
where SDF (in %) is the stream depletion factor of the model cell to which the new well is added; Qr0 and Qr are the exchange rates between streams and aquifers after and before adding the new pumping well in the model; and Q p is the pumping rate of the added hypothetical well. As shown in the equation, the streamflow depletion factor was estimated through modeling the effect of a new hypothetical well pumping continuously at a constant rate for 50 years based on the baseline model. In this study, the calibrated LPR groundwater model serves as the baseline model. Conventionally, the depletion estimation of each hypothetical well needs an individual simulation by adding it to the baseline model. The hypothetical well initially pumps the water from groundwater storage. As pumping continues, the drawdown caused by the new pumpage increases and expands, thus increasing the outflow from the stream or reducing the inflow to the stream when the induced drawdown reaches the stream. Eventually, the depletion rate becomes stable as the drawdown
55 reaches other water sources such as other streams, general head boundaries or constant head boundaries and the system reach another steady state. To account for the depletion change with time, it needs to run both the baseline model and the depletion model in the transient state. Though the groundwater pumping rate can quantitatively change the exchange rate between streams and aquifers, it have proved that the pumping rate of the hypothetical well has a relatively small effect on the stream depletion factor. However, large pumping rate will dewater the grid cell rapidly. Because the estimation of stream depletion needs to be performed at a large number of model grid cell, it can be very computational costly to carry out such analysis. Therefore, a simpler but effective stream depletion program Stream Depletion Analysis (SDA) is developed to linearize the flow equations.
2.6.1
Explanation of MODFLOW-SDA
In the analytical solution developed by Hunt (1999) for stream depletion, SDF is irrelevant to the pumping rate (Q p ) used for the new pumping well. The volume in each source contributing to the pumping will change with Q p . However, the proportional distribution only changes with the location of the well, which is attributed to the ideal assumptions used in the solution. In Hunt’s solution, it’s assumed that (1) drawdowns are small enough compares with the saturated aquifer thickness to allow the governing equations to be linearized; (2) Changes in water surface elevation in the river created by pumping are small compared with changes created in the water table elevation on the aquifer side of the semipervious layer. Analogically, such assumptions are made in MODFLOW to linearize the flow equation. In MODFLOW, there are generally two types of flow boundary
56 conditions according to the flow’s dependency on groundwater heads. The flows of head-independent boundary conditions, which include wells and recharges, are specified directly as model inputs. In contrast, the flows of head-dependent boundary conditions need to be calculated after the groundwater flow equation is solved. The general equation of the head-dependent boundary conditions can be expressed as:
Q HB = CB ( HB − h)
(2.13)
where Q HB is the flow from the boundary to the aquifer; CB is the flow conductance coefficient; HB is the boundary head which, for example, may be the stream stage or ET extinction depth; and h is the groundwater hydraulic head. MODFLOW employs a finite-difference method to calculate groundwater heads and flows. Converting the groundwater flow equation into finite-difference form yields: CCi−1/2,j,k hi−1,j,k + CRi,j−1/2,k hi,j−1,k + CVi,j,k−1/2 hi,j,k−1
+(−CCi−1/2,j,k − CRi,j−1/2,k − CVi,j,k−1/2 − CCi+1/2,j,k − CRi,j+1/2,k − CVi,j,k+1/2 + HCOFi,j,k )hi,j,k +CCi+1/2,j,k hi+1,j,k + CRi,j+1/2,k hi,j+1,k + CVi,j,k+1/2 hi,j,k+1 = RHSi,j,k (2.14)
where CC, CR and CV are the conductance along the X, Y and Z directions; HCOF and RHS are coefficients related to the sources and/or sinks and storage terms; h is the hydraulic head; the subscripts i, j and k are used to designate the cell column, row and layer respectively; and 1/2 denotes the region between two grid cells. The boundary stresses are added to or subtracted from the groundwater
57 storage through the HCOF and RHS terms: N1
HCOF = −
S
∑ CBn − ∆tC
n =1 N2
RHS = −
∑
n =1
QB n −
N1 SC h0 − ∑ CB n HB n ∆t n =1
(2.15) (2.16)
where QB is the head-independent boundary flow; SC is the storage capacity; ∆t is the time step length; N1 and N2 are the numbers of head-dependent and head-independent boundaries. In a scenario run, the change of the flow system is usually justified by changing the RHS term. For example, adding pumping wells or reducing groundwater recharge lead to a reduced QB value. The entire system of the flow equations can be expressed in the matrix form as:
Ah = q
(2.17)
A0 h0 = q0
(2.18)
where A is a matrix of coefficients of heads; h is the groundwater head at each cell; q is the initial storage and direct flow terms such as recharge and pumping; the superscript prime indicates the scenario run. Matrix A and A0 are related to all the head coefficients and boundary conductances based on Equation 2.14 and Equation 2.15. For each solver iteration A is re-formulated based on the updated head values. SC does not change with heads unless the cell is converted between unconfined and confined conditions. In the scenario run, therefore, the change of A includes deviations of CC, CR, CV and CB , which can be attributed to the new pumping stress. CC, CR, CV depend only on layer thickness in confined layers. For an unconfined layer, these coefficients become head dependent, resulting in
58 changed values in the scenario run. In a regional groundwater model, however, the head changes between two runs are usually minimal compared with the aquifer thickness. Therefore, the first assumption of the analytical solution can also be applied to MODFLOW. The changes of CC, CR and CV, caused by the head changes, can be neglected. Meanwhile in the stream package, the streamflow will change as the baseflow has been altered due to the groundwater head differences. Accordingly, CB can be changed as the wet radius changes with the stream stage which is related to the streamflow. With the assumption that the accretion flow is negligible to the streamflow, the change of the stream stage can also be ignored. Under the same assumptions (1) and (2) of the Hunts solution, the change of A due to the groundwater head change is neglected. Therefore, the flow equation is linearized and we can subtract Equation 2.18 from Equation 2.17: A( h0 − h) = (q0 − q)
(2.19)
A∆h = ∆q
(2.20)
which yields:
where ∆q = −∆Q B − ∆S0 . ∆Q B represents the specified accretion flow. The change of initial storage ∆S0 is induced by the head change in the previous time step (∆h0 ): ∆S0 =
Sc ∆h0 ∆t
(2.21)
Clearly in the first time step, ∆S0 = 0 as the initial heads between the two model runs are the same. After linearization of A, it becomes self evident that the change of head is proportional to the accretion flow and therefore SDF is not related to the accretion flow rate. This, in turn, explains that how SDA results can be effected by
59 different pumping rates of the hypothetical well. Equation 2.21 can be solved by using one of the MODFLOW solvers. Since A has been calculated in the baseline model and it does not change between solver iterations, the solution becomes very efficient and it doesn’t require to read the MODFLOW input files again in the scenario run. The change in the head-dependent boundary flow and storage can therefore be calculated using ∆h: ∆Q HB = −CB ∆h ∆S = −Sc (∆h − ∆h0 )
(2.22) (2.23)
Based on the linearized equation, a package is developed for stream depletion analysis in MODFLOW (MODFLOW-SDA). The flow chat of the package is shown in Figure 2.27. In the baseline run, the coefficient values are extracted are stored for the scenario runs. IBOUND, which represents the cell type, is needed to define the constant head boundary. The type, location and contuctance of each head-dependent boundary are also saved for calculating the flow difference based on the head change. MODFLOW-SDA is not restricted to streamflow depletion calculation. It also calculated the flow difference of other stress types. To improve the read/write efficiency and reduce storage space, the extracted information from the baseline run are saved as binary files.
2.6.2
Application of MODFLOW-SDA
MODFLOW-SDA is applied to a simple hypothetical model to assess its accuracy. After the verification on the hypothetical model, it is used on the LPR model.
60
Figure 2.27: Flow chart of MODFLOW-SDA. 2.6.2.1
Hypothetical case
The grid of the hypothetical model is composed of seven columns, seven rows and one layer with uniform a cell size. The row and column widths of a cell are 100 m and the layer thickness is 50 m. The River, Constant Head and Well packages are used to represent the boundary conditions. The river is located at the top row which consists seven cells. There is one constant head cell at the lower left conner. The pumping well is located at the center of the model which is the fourth column and row. At the beginning of the simulation, the river stage, groundwater initial head and the constant head are at the same level, representing a system of equilibrium (no flow condition). The well starts to pump in the first stress period which lasts for 10 days and is deactivated in the second period for another 10 days. The stream leakage after pumping starts is considered as the depletion compared
61 to the equilibrium condition. In this test, different pumping rates are used and the results of the SDA method (SDA) and the two run method (TRM) are shown in Figure 2.28. It’s easy to observe from Figure 2.28 that the SDF does not change with pumping rates unless the pumping rate is too large and convert the well cell to a dry cell. This occurs when pumping rate is 105 m3 /d. Meanwhile, the head changes (drawdown) have very little effect on SDF results. The maximum difference of SDFs between SDA and TRM occurs at the end of the pumping with the rate of 10000 m3 /d. The drawdown reaches 8.15 m which is over 15 % of the thickness, indicating that the transmissivity is reduced over 15 %. The error in the SDA model, however, is only 0.3 %. 2.6.2.2
LPR model case
The groundwater flow model can predict changes of the exchange rate between streams and aquifers. Through adding new wells at different locations in the modeling area, the SDF of a specific cell can be calculated. The stream is consider to be hydrologically connected to the aquifer if the depletion ratio of the new well maintains lower or equal 10% at the end of a 50-year continuous pumping period. The hypothetical wells are set up according to the following steps: 1. Hypothetical wells are added along the Platte River and subsequently outwards until SDF was reduced to lower 10%. 2. Hypothetical wells are screened in the first layer and operated at a constant pumping rate of 100 GPM. A total of 122 hypothetical wells were added to the model for the streamflow depletion analysis. Table 2.5 shows the SDF for selected hypothetical wells. The
62
Pumping rate: 1 m3/d; Drawdown: 0.000757938 m Pumping rate: 10 m3/d; Drawdown: 0.00757961 m 15
●
●
● ●
●
● ●
●
●
●
10 ●
●
5 ●
●
●
●
0
Stream depletion factor, %
Pumping rate: 100 m3/d; Drawdown: 0.0758395 m Pumping rate: 1000 m3/d; Drawdown: 0.762778 m 15
●
●
● ●
●
● ●
●
●
●
10
Model ●
●
●
SDA TRM
5 ●
●
●
●
0 Pumping rate: 10000 m3/d; Drawdown: 8.14758 m Pumping rate: 1e+05 m3/d; Drawdown: −1e+20 m 15
●
●
● ●
●
● ●
●
●
●
10 ●
●
5 ●
●
●
●
0 5
10
15
20
5
10
15
20
Time
Figure 2.28: SDF results of the hypothetical model. individual hydrographs of the simulated SDFs with two methods are shown in Appendix A. Hereafter, the error refers to the difference in the calculated SDFs between SDA and TRM.
63 Table 2.5: Locations and depletion rates of hypothetical wells. Row
Column
SDF
Mean error
18 19 20 20 20 20 20 22 22 22 23 23 25 25 25 25 25 25 25 25 25 26 26 27 27 27 28 28 28 30 30 30 30 30 30 30 30
95 95 75 85 95 105 110 60 85 110 60 110 60 75 85 90 95 105 120 130 140 120 130 120 130 140 120 130 140 50 60 90 105 110 120 130 140
23.25% 26.07% 27.20% 27.77% 28.52% 18.93% 12.29% 0.64% 34.41% 13.51% 19.93% 13.68% 20.51% 37.51% 59.33% 63.49% 52.24% 34.21% 0.20% 0.28% 0.24% 7.00% 5.37% 22.38% 10.44% 3.58% 35.91% 15.90% 7.33% 7.26% 25.13% 88.23% 60.85% 58.59% 58.14% 25.59% 15.60%
4.37% 3.98% 6.90% 7.68% 3.64% 0.55% -0.51% 17.32% 7.86% 0.18% 5.32% 0.47% 4.48% 5.47% 5.31% 3.46% 2.87% 0.99% -0.01% -0.02% -0.02% -0.58% -0.47% -1.29% -0.98% -0.17% -1.62% -1.59% -0.36% 2.30% 2.29% 1.26% -0.32% -0.42% -1.89% -2.62% -0.89%
Continued on next page
64 Table 2.5 – continued from previous page Row
Column
SDF
Mean error
35 45 48 49 49 50 50 50 50 50 51 52 53 55 56 63 64 65 70 70 98 98 99 100 100 100 100 100 100 100 100 100 101 101 103 103 103 104
50 75 75 10 75 20 30 40 75 85 10 172 40 30 170 3 3 30 85 130 105 140 105 3 10 30 75 90 95 105 130 140 10 144 3 10 30 3
7.94% 42.07% 19.54% 4.62% 11.78% 7.91% 6.27% 5.97% 2.30% 82.08% 7.05% 9.47% 7.80% 5.34% 15.72% 7.37% 8.22% 38.54% 31.00% 84.58% 12.52% 5.46% 5.41% 33.24% 10.73% 14.19% 3.04% 6.28% 11.53% 1.84% 2.20% 2.69% 9.68% 11.43% 15.56% 7.56% 11.80% 11.57%
1.11% -3.27% -2.83% -0.08% -1.95% 1.16% -0.30% 0.06% -0.40% 2.14% -0.14% -0.80% 0.03% -0.34% -1.39% 1.05% 1.19% -2.04% 1.49% 3.82% 1.09% 0.93% 0.44% 1.55% 0.36% 3.50% 0.34% 0.52% 2.41% 0.14% 0.29% 0.58% 0.42% 1.65% 1.10% 0.49% 3.21% 0.94%
Continued on next page
65 Table 2.5 – continued from previous page Row
Column
SDF
Mean error
104 104 105 105 105 105 105 105 106 107 107 107 110 110 110 110 110 111 111 112 112 114 115 115 116 117 124 125 126 127 140 140 140 150 150 150 150 150
30 142 3 10 30 50 95 141 95 3 50 95 3 30 95 130 140 198 200 140 196 143 40 140 139 138 144 140 139 138 155 158 160 160 165 166 167 170
10.35% 13.16% 8.00% 6.28% 9.15% 7.89% 2.77% 12.22% 2.05% 2.04% 7.47% 1.49% 1.27% 4.09% 0.51% 2.33% 6.72% 10.78% 3.65% 4.29% 18.10% 18.32% 12.72% 4.48% 4.44% 4.43% 14.86% 7.04% 5.39% 3.70% 4.88% 10.41% 15.50% 3.62% 8.24% 11.90% 15.75% 27.84%
2.91% 2.23% 0.75% 0.52% 2.65% -1.25% 0.43% 2.01% 0.27% 0.27% -1.27% 0.17% 0.25% 1.38% 0.04% 0.38% 1.08% -0.97% -0.32% 0.77% -1.52% 2.60% 4.71% 0.81% 0.81% 0.82% 3.16% 1.60% 1.24% 0.87% 0.11% 0.08% -0.02% -0.01% -0.18% -0.24% -0.30% -0.50%
Continued on next page
66 Table 2.5 – continued from previous page Row
Column
SDF
Mean error
160 160 160 160 160 170 170 170 179 179
160 165 168 170 180 170 175 180 194 195
3.00% 3.43% 11.46% 18.61% 60.33% 2.88% 12.93% 32.76% 12.40% 20.42%
0.08% -0.02% -0.20% -0.39% -0.87% 0.01% -0.01% -0.25% 0.08% -0.15%
As shown in the distribution of mean errors (Figure 2.29), overall, the results of SDA are favorable with those calculated with TRM. The mean errors are less than 1% for over 50 % of the wells. For large discrepancy, they are caused by the numerical noises in TRM. In the low permeable area, where pumping can create large drawdowns, the errors are also larger.
2.6.3
Discussion on MODFLOW-SDA
There are some benefits by using SDA over TRM: 1. SDA can avoid dry cells when the water supply capacity of the cell is lower than the pumping rate of the hypothetical well. 2. The computer time can be dramatically reduced because of the more efficient algorithm and input/output processes. For the LPR model, an average of four fold speedup has been achieved. 3. The head change becomes the solver target in SDA, which to some extent lessen the numerical noises inherent in TRM.
67
density
0.4
0.2
0.0 0.0
2.5
5.0
7.5
Mean error, % Figure 2.29: Distribution of mean errors of SDF results simulated with SDA.
2.7 Summary The regional ground-water flow system in the Lower Platte River basin was simulated by a three-dimensional model for the time period from 1950 to 2004. The model was constructed with MODFLOW-2005, a three dimensional, finitedifference, ground-water flow model. The model was calibrated through matching the observed data of 189 monitoring wells. The model achieved high accuracy of simulating hydraulic heads, therefore providing a reliable tool for the regional groundwater flow and streamflow depletion analysis. Groundwater recharge was dominated by precipitation and irrigation return flow. The modeled recharge was estimated through the infiltration ratio. The groundwater pumpage and river leakage composed the most discharge from the
68 groundwater system. Other terms such as groundwater evapotranspiration and boundary flux were also simulated but they were relatively low compared to the main terms in the water budget of the groundwater model. Streambed data were collected in the field. The stream width and depth were estimated through empirical functions. The calibrated groundwater flow model was used to predict the irrigation effects on the streamflow of the Lower Platte River. The exchange rate between streams and aquifers were estimated based on the difference of stream stage and hydraulic head of the underlying aquifer. Increasing pumpage can induce the change of stream-aquifer exchange rates through lowering the hydraulic heads. This effect is more significant as the well is closer to the stream. During the modeling period from 1950 to 2004, the groundwater level did not decline much in the study area and streamflow records did not show a decline trend for most gauging stations in the study area (Wen and Chen, 2006). A stream depletion analysis tool has been developed to evaluate the stream-aquifer connectedness. SDA is effective on stream depletion analysis and reduces numerical noises. It also eliminates the effect of pumping rate used for the hypothetical wells. The model documented in this report has some limitations. Uncertainties in some model inputs are difficult to quantify. Groundwater recharge, for example, is comprised by precipitation and irrigation return flow, but the infiltration processes are substantially complicated. Groundwater evapotranspiration is also a complicated component in the hydrologic cycle. Though, the model still achieved reasonable accuracy compared to the measurements of hydraulic heads, integrated modeling tool can be used to represent the complex processes of surface water and groundwater interaction.
69
Chapter 3 A CROSS-SECTION BASED STREAMFLOW ROUTING PACKAGE FOR MODFLOW 3.1 Introduction Understanding hydrological processes that occur in the stream-aquifer systems has significant importance to integrated water management, especially where streams and aquifers are hydrologically connected (Rassam, 2011; Sophocleous, 2002; Winter, 1995). The dynamics of stream-aquifer water exchanges are also important to the functioning of aquatic ecosystems (Hancock et al., 2005), pollutant and nutrient transport (Findlay, 1995; Jones and Holmes, 1996), as well as the quality and quantity of water supply to domestic, agricultural and recreational purposes (Winter, 1995). MODFLOW is the most widely used numerical model to simulate stream and groundwater interactions. However, it is challenging to represent stream-aquifer interactions for wide channel streams or for high-resolution grids, especially
70 when the heterogeneity of streambed properties across the channel and stream morphology need to be taken into account. The spatial variability of water and solute fluxes at the groundwater-surface water interface, the hyporheic zone, and their impacts on the stream and riparian ecosystems were studied extensively through experimental measurements in previous studies and continue as an area of active investigation (Biksey and Gross, 2001; Binley et al., 2013; Faulkner et al., 2012; Rosenberry and Pitlick, 2009; Wondzell et al., 2009). Such variability is largely attributed to streambed heterogeneity which is universal because of the stream flow dynamics and movement of streambed sediments, as well as streambed chemical reactions and invertebrate activities (Dong et al., 2012). Chen et al. (2008) discovered high spatial heterogeneity in streambed materials with depth, along the channel and cross the channel in the Platte River, Nebraska. Lu et al. (2011) also demonstrated that distribution of both horizontal and vertical hydraulic conductivities is significantly spatially dependent. Heterogeneity in permeability of streambed materials can form many flowpath connections between streams and aquifers, for both water and solutes (Conant et al., 2004) , and consequently provide habitats of various scales which are important to different species (Brunke and Gonser, 1997). It has been demonstrated that neglecting streambed heterogeneity could cause significant errors in modeling results of stream-aquifer interactions (Irvine et al., 2012). Brunner et al. (2010) pointed out the limitation in horizontal discretization of rivers could result in biased estimation of stream-aquifer exchange flux when the river width and cell size are mismatched. Mehl and Hill (2010) also examined the effect of the MODFLOW grid size on the simulation of stream-aquifer interactions. In the streamflow depletion analysis. On the other hand, Chen et al. (2008) showed that the depletion rate induced by groundwater pumping can be affected by
71 the grid cell sizes. They all proved that refining the model grid can improve simulation accuracy. These findings underline the need for models to represent the geometry and properties both across and along stream channels with grids of high resolutions. Because of the importance of streambed heterogeneity, and, particularly the need for simulating stream-aquifer interactions in wide channels or with refined grids, the Cross-Section based streamflow Routing (CSR) package was developed for MODFLOW. The objectives for developing the CSR package were: (1) to simulate stream-aquifer interaction and streamflow routing for the wide channel streams or with high-resolution grids; (2) to represent the heterogeneity in streambed properties not only along stream channels but also across them; (3) to physically represent the streambed geometry of high complexity. In this paper, we describe the concepts and algorithms used in the CSR package. A hypothetical stream-aquifer system is used to test the capabilities of the package. We also discuss the capabilities, assumptions and limitations, as well as potential improvement of the package.
3.2 Literature review on River and Stream Packages of MODFLOW MODFLOW (McDonald and Harbaugh, 1988) was written in FORTRAN 77 and divided specific hydrologic features and calculation functions into separate modules. The modular structure, providing ease of modifying and adding features, is still being used in the current version. MODFLOW numerically solves the following finite-difference equation that describes the three-dimensional groundwater flow
72 (Equation 2.14). A River (RIV) package was incorporated in the first version of MODFLOW. According to Darcy’s law, Walton and Ackroyd (1966) and Prickett and Lonnquist (1971) indicated that the seepage between streams and aquifers could be calculated based on the streambed area, head difference between streams and aquifers and the streambed permeability. Similar to their research, the RIV package assumed that the head loss between surface water and groundwater occurred across their interface (Rushton, 2007), which is commonly considered as the streambed. The conductance (CDT), a streambed coefficient representing streambed’s capacity of transmitting water, was incorporated in the following equation to calculate the seepage Qgwc of a grid cell:
Qgwc =
CDT( Hs − Hg ),
Hg ≥ BOT
CDT( Hs − BOT),
Hg < BOT
(3.1)
where Hs and Hg are the water levels of the stream and groundwater, respectively; BOT is the bottom elevation of the streambed; CDT = Kb W L/M; Kb is the streambed hydraulic conductivity; W and L are the stream width and length in the grid cell, respectively; and M is the streambed thickness. This method has been widely applied to calculate the seepage to or from surface water in other MODFLOW packages. In the RIV package, the river stages are specified by users and remain constant during a model stress period, which implies that the river can act as an infinite source or sink. Nevertheless, this assumption may flatly contradict the reality under some circumstances. The RIV package is not able to correctly represent the stream-aquifer interaction where the stream flows and stages, for example in headwater sections that are fed by groundwater baseflow,
73 are sensitive to the streambed seepage. Furthermore, concerns of the influence of streamflow depletion induced by pumping wells had increased. There was a need for a model that can deliver accurate representations of both stream-aquifer interactions and water flow in streams. Subsequently, Prudic (1989) developed the Stream (STR) package which was a tremendous improvement compared to the RIV package. The most profound change of STR is that it could simulate the streamflow routing in complicated stream networks which allow diversions and conjunctions. In the STR package, the streambed seepage calculation is also performed based on Equation 3.1. However, seepage to aquifers is restricted when the stream reach becomes dry. The stream stage is solved based on the relationship of the streamflow and water depth described by Manning’s equation. Rivers and streams are both represented by rectangular channels in the RIV and STR packages. To simulate contaminant transport within a stream network and between streams and aquifers, a new Streamflow-Routing (SFR1) package was developed to improve the STR package (Prudic et al., 2004). The package can be linked with the Lake (LAK3) package (Merritt and Konikow, 2000) and the Ground-Water Transport (GWT) Process package (Konikow et al., 1996). It allows users to specify the overland runoff, precipitation, evaporation at each stream cell. The SFR1 package also provides additional options to calculate the stream water depth and representation of eight-point cross-sections. The SFR2 package, the successor of SFR1, restricts the unsaturated flow under the streambed when the groundwater level drops below it (Niswonger and Prudic, 2005). This function was facilitated by simulating the unsaturated flow beneath streams using the same method as the UZF package (Niswonger et al., 2006). The streamflow routing method in SFR1 and early versions of SFR2 is identical to the STR package by using the continuity equation without accounting for flow attenuation. However, SFR2 was enhanced
74 in GSFLOW (Markstrom et al., 2008) through adding the kinematic wave routing approach. Officially incorporated in MODFLOW, the RIV, STR and SFR packages are the most widely used packages to simulate stream-aquifer interactions. With the development of computer programs in stream routing models, on the other hand, more complicated stream routing models were directly coupled with MODFLOW. MODBRANCH (Swain, 1993), combining MODFLOW and a more sophisticated stream routing model BRANCH (Schaffranek et al., 1981), was released to simulate unsteady flow in open-channel networks. By solving the St. Venant equation, MODBRANCH can overcome the limitation of the STR package that cannot simulate the transient and non-uniform flow conditions. Accounting for difference in the time scales of the surface water flow and groundwater flow, MODBRANCH also allows streamflow routing to be computed in shorter lengths of time intervals than the groundwater flow computation. DAFLOW-MODFLOW, coupling DAFLOW (Jobson, 1989) and MODFLOW, is also able to simulate the transient flood movements. Differing from BRANCH, DAFLOW solves a simplified form of the dynamic wave equations by assuming unidirectional flood wave movement. It is therefore, not able to simulate the backwater or bidirectional flow conditions. According to the authors of DAFLOW, however, it could provide a relatively stable solution that requires minimum field data and calibration. DAFLOW was incorporated in MODFLOW-2000 as the DAF1 Package (Jobson and Harbaugh, 1999). The Surface-Water Routing Process (SWR1) package was recently developed to simulate one-dimensional and two-dimensional surface-water flow (Hughes et al., 2012) for MODFLOW-2005. SWR1 implements the reservoir routing and diffusive-wave flow routing methods. It can account for backwater (tailwater) effects, bidirectional surface-water flow, and management of surface water using control structures.
75 Other attempts have been made by modifying the river or stream packages in MODFLOW or coupling other streamflow models with MODFLOW to meet particular modeling needs. MOBFLOW, developed by Osman and Bruen (2002), can account for the effect of the unsaturated flow conditions where the groundwater levels fall below streambeds by modifying the equations of the RIV package. Rodriguez et al. (2008) combined HEC-RAS (Brunner, 1995) to simulate the interaction between the drainage network and aquifers. To evaluate the effect of river and reservoir operations, Valerio et al. (2010) linked MODFLOW with the object-oriented river and reservoir operation model, RiverWare. In addition, other integrated surface water and groundwater flow models such as SFWMM-MODFLOW (Yan and Smith, 1994), FHM (Ross et al., 1997), SWATMOD (Sophocleous et al., 1999; Kim et al., 2008) and GSFLOW (Markstrom et al., 2008) coupled the watershed hydrological models with MODFLOW to simulate the dynamics of surface water, groundwater, and stream-aquifer interactions. Table 3.1 compares the capabilities between nine stream packages developed for MODFLOW. The streams are represented as rectangular channels in nearly all the packages except for SFR and SWR1 in which cross-sections can also be represented by eight streambed points or irregular geometries, respectively. RIV, RiverWare and SWR1 are the only models capable of representing the channel width spanning over multiple cells. Nevertheless, the RIV package does not provide flow routing simulation. RiverWare cannot represent complex cross-section geometries. SWR1 requires grouping a collection of model cells in two-dimensional geometric forms to physically represent wide channels spanning over multiple cells, in which case the geometry of the cross-section is determined by the top elevation of each cell. In other packages, wide channels can also be represented by multiple parallel stream segments but this method increases the modeling work and may result in
76 Table 3.1: Comparison of MODFLOW stream packages. Package
Routing Method
Geometry
Wide Channel
Unsaturated Flow
RIV
Null
Rectangular
Y
N
STR
Mass Conservation
Rectangular
N
N
SFR1
Mass Conservation
Rectangular Eight-Point
N
N
SFR2
Mass Conservation Kinematic Wave
Rectangular Eight-Point
N
Y
SWR1
Reservoir-Routing Diffusive Wave
Rectangular Trapezoidal Irregular
Y
N
BRANCH
St. Venant’s Eq.
Rectangular
N
N
DAF1
Diffusive wave
Rectangular
N
N
HEC-RAS∗
St. Venant’s Eq.
Rectangular
N
N
RiverWare
Time lag Rectangular Y N Muskingum-Cunge Kinematic Wave ∗ HEC-RAS can account for irregular cross-section geometries. The stream (drain) cross-sections were assumed as rectangular in Rodriguez et al. (2008). uneven stream stages between the segments at the same cross-section. Besides, it is difficult to represent the variation of streambed hydraulic properties across the cross-section. Development of the CSR package is expected to overcome these limitations and thus improve the simulation accuracy of stream-aquifer interactions.
77
3.3 Development of the CSR package 3.3.1
Overview of the model
The CSR package is developed to incorporate the capability of simulating streamflows and stream-aquifer interactions where the stream width is larger than the MODFLOW grid cell size. In the CSR package, the stream components are represented by different conceptualized geometry elements. Instead of a segment of polyline employed in most other stream packages to represent a stream segment, CSR uses a four-point polygon (quadrilateral) that is composed of an upstream cross-section and a downstream cross-section. A cross-section consists of a number of streambed points possessing properties that describe the streambed geometry and hydraulic properties. The left and right end points are used to determine the locations of the stream segments. Based on the cross-section geometry and hydraulic properties, CSR calculates the new stream stage at the downstream cross-section using Brent’s method to solve Manning’s Equation. A module is developed to automatically compute the submerged area of the stream segment polygon on each intersected MODFLOW grid cell as the upstream and downstream stages change. The stream stage and streambed hydraulic properties of a model cell are interpolated based on upstream and downstream stages and the streambed points, respectively. Streambed leakage is computed as a function of the streambed conductance and difference between the groundwater level and stream stage. The mass conservation method or the Muskingum-Cunge flow routing scheme with variable parameters are used to simulate the streamflow as the groundwater (discharge or recharge) contributes as lateral flows.
78
3.3.2
Conceptualization of stream components
The CSR package has four types of conceptualized elements representing stream components, namely streambed points (SBP), cross-sections (CS), stream segments (SEG) and stream networks (NET). Each SBP carries the spatial information and hydraulic properties that are needed in streamflow calculation. A CS is formed by an arbitrary number of SBPs that are capable of representing the geometry and the variability in hydraulic properties of the CS (Figure 3.1). The example CS shown in Figure 3.1 is defined by five SBPs. The SBPs at the two ends of the CS are the full-bank points (FBPs) that are always higher than the stream stage during the whole simulation. As shown in Figure 3.1, SBP [a] and [e], for example, act as the FBPs. To determine the horizontal position of the CS in the model, users need to define the model coordinates (x, y, z) of the FBPs. Other SBPs such as [b], [c] and [d] only need to be specified with their horizontal and vertical distance (HD, VD) to the left FBP. As indicated in Figure 3.1, the location of stream-aquifer interaction varies with fluctuation of the stream water surface. For example, the submerged streambeds span over the Column 2 – 9, 2 – 8 and only 6 – 7 when the river stages are at the positions of h1, h2 and h3, respectively. Two adjacent CSs, an upstream CS (UPCS) and a downstream CS (DNCS), constitute a SEG. In other words, SEGs are connected by a shared CS. For instance CS [b-b’] in Figure 3.2, is the DNCS of SEG [I] whereas it is also the UPCS of SEG [II]. Besides the shared CSs, there are also headwater CSs [a-a’] and [d-d’] and outlet CSs [c-c’] and [e-e’] (Figure 3.2). Users need to specify the inflow of the headwater CSs at each stress period. Tributary SEGs could join to the main channel SEGs by adding their outlet streamflows as lateral flows to the main channels. The SEGs must be numbered sequentially from upstream to downstream because
79 f
g
a
e
h1 c
1 h2
2 3
h3 b
4
d
stream cross-sectoin streambed point MODFLOW grid cell
5 1
2
3
4
5
6
7
8
9
10
Figure 3.1: Schematic diagram of an example MODFLOW grid superimposed by a stream segment. the SEG number determines the order of streamflow calculation. The number of the tributary SEG also needs to be higher than the main channel SEG to which it connects. Hence, the outflows from upstream SEGs can be computed before starting the streamflow calculation of the downstream SEG and will be correctly added as inflows. In contrast, the numbers of CSs are not important as long as the corresponding numbers of the UPCS and DNCS are properly assigned to each SEG. Diverting flow can also be simulated with diverting SEGs. The diverted flow is restricted not to exceed the upstream flow of the SEG from which the water is diverted.
3.3.3
Computation of stream depth
The stream depth during each model stress period can be directly assigned by users or computed based on the relationship with streamflow. There are four methods in CSR to define the relationship between the stream depth and streamflow. The first two options are both based on Manning’s Equation which is widely used in
80 a
d
d’
a’
II b
I
e
e’ c
b’ III c’ full-bank point
I, II and III stream segment number
submerged area
stream segment
Figure 3.2: Plan view showing an example MODFLOW grid superimposed by three connected stream segments. hydraulic engineering and surface hydrology: C Qs = ( ) AR2/3 S1/2 n
(3.2)
where Qs is the stream volumetric flow rate; C is the unit factor, which is equal to 1 or 1.486 for metric units or imperial units, respectively; A is the cross-sectional flow area below the stream stage; n is Manning’s roughness coefficient which could be estimated based on reference values (Chow, 1959); R is the hydraulic radius of the stream and S is the streambed slope. As the stream stage changes, the flow area will vary accordingly and it is calculated based on the stream stage and streambed geometry (locations of SBPs). Due to the variability in the roughness and slope between SBPs, a CS can be separated into subdivisions according to the SBP locations (Figure 3.3). The representative values of n and S are computed through linear interpolation at the midpoint of the submerged streambed section
81
stream stage I
II
III
IV stream depth
midpoint location streambed
Figure 3.3: Cross-section subdivision for stream depth calculation. in each subdivision. The streamflow is subsequently computed through the Single Channel Method (SCM) or Divided Channel Method (DCM) in Option 1 and 2, respectively. In the SCM, the representative values of n and S are subsequently averaged by weighting them according to their subdivision wet perimeters:
nc =
∑(ni Pi )/ ∑( Pi )
(3.3)
Sc =
∑(Si Pi )/ ∑( Pi )
(3.4)
where Pi , ni and Si represent the wet perimeter, representative n and S of the subdivisions; and the subscript i is the subdivision index. After computing the composite values nc and Sc , they are substituted into Equation 3.2 to achieve the CS streamflow. Alternatively, the DCM separately calculates the sub-flow of each subdivision using Manning’s equation and subsequently sums up the sub-flows:
Qs =
1
1/2 ∑(Qsubi ) = C ∑( ni Ai R2/3 i Si )
(3.5)
82 where Qsubi is the flow of the subdivision; Ai is the subdivision area and Ri = Ai /Pi . If transverse variations in n and S are insignificant, the SCM performs better than DCM and vice versa. Because of the complexity of the relationship between stream depth and streamflow in the compound channels, the stream depth cannot be explicitly solved from the equation. An improved Brent’s method (Zhang, 2011) is applied to solve the equations numerically. The third option uses a power function to approximate the change of the stream depth with the streamflow:
Qs = αDs β
(3.6)
where Ds is the stream water depth; and α and β are constant coefficients specified by users. The fourth option allows users to design the relationship between the stream depth and streamflow by entering data entries of streamflows and the corresponding water depths. CSR will linearly interpolate the water depth based on the two closest data entries embracing the desired streamflow. If the streamflow is larger than the maximum tabulated flow rate, the depth is extrapolated from the two maximum data entries. After the computation of stream depth, CSR computes average flow velocity and the width of the submerged streambed section. Consequently, it can determine locations of the intersected model grid cells where the stream-aquifer interactions actually occur.
83
3.3.4
Computation of intersections of the submerged area
Based on the upstream and downstream submerged streambed locations, a routine is incorporated into CSR to rapidly compute the intersected submerged area on each single grid cell. The routine is a modification based on the method developed by Larkin (1988). An example is demonstrated in Figure 3.4 to explain the algorithm of this method (referred as SBAREA hereafter) to compute the intersected area. The algorithm requires that the vertices of the submerged quadrilateral are arranged in the clockwise orientation. This premise is ascertained by the SORTCLOCK4 function of the CSR code. 1
1
2
2
3
3
4
4
5
5 1
2
3
4
5
6
Step 1
1
2
3
4
5
6
Step 2 Symbol explanation
1
searching point
2
submerged area
3
added area
4
subtracted area
5 1
2
3
4
5
6
overlapping area
Step 3 Figure 3.4: Algorithm of computation of the intersection between the submerged streambed area and the MODFLOW grid (plan view).
In Step 1, SBAREA searches for the intersected points between the edges of
84 the quadrilateral and the right vertical side of each grid cells. If the edge passes through the cell side from left to right, the cell is along the top of the quadrilateral, and the area of the rectangle under the intersected point (hatched as diamonds in Figure 3.4) will be added to the cell area. Contrarily, the cell area will be subtracted by the area of the rectangle under the intersected point (hatched as diagonal lines in Figure 3.4) if the edge penetrates from right to left. In this case, the cell is cut by the bottom edge of the submerged quadrilateral. Step 2 is to search for the cells with their sides passed through by the edges of the quadrilateral. Two intersected points form a trapezoid or triangle with the left vertical side of the cell. If the edge starts or ends within the cell, its end point acts as another intersected point. The areas of the trapezoids and triangles are added to or subtracted from the cell areas when the edges penetrate the cell downwards or upwards, respectively. In Figure 3.4, the overlapping area, shaded in black, represents a part of the area that is covered by both the added and subtracted areas, for example, at the top left cell and the bottom right cell. The third step is to include the whole area of the cell whose upper right corner is contained by the quadrilateral. These cells can easily be located through the Ray Casting Algorithm which is the simplest method of finding whether a point is inside or outside a polygon. The algorithm counts the number of intersections between the polygon and a ray starting from the point. The point is inside the polygon if the number of intersections is an odd number, and it is outside if it is an even number. The final result is the summation of the subtracted and added areas of each MODFLOW grid cell. In fact, all the steps are executed simultaneously on each grid cell when it searches along the edges of the quadrilateral. Therefore, the intersected submerged area on each MODFLOW grid cell can be efficiently computed with
85 SBAREA by scanning the perimeter of the submerged area only once.
3.3.5
Computation of the stream-aquifer seepage
The method used to compute the stream-aquifer flux is similar to that used in other stream packages of MODFLOW. However, instead of using the width and length of the stream reach, the actual submerged area is calculated prior to the computation of stream-aquifer seepage and is used to calculate the streambed conductance:
CDT = Kb A0 /M
(3.7)
where A0 is the submerged area projected on the horizontal plane of a MODFLOW grid cell. The stream-aquifer seepage in a cell is consequently calculated using Equation 3.1 with the computed value of CDT. This method neglects the horizontal components of the stream-aquifer seepage occurring on the sloping streambeds. The occurring layer of the interaction is determined based on elevation of the streambed bottom and the vertical discretization of the grid cell. The stream stage at the grid cell is calculated through linear interpolation of the stream stages at the UPCS and DNCS:
d=
|( x2 − x1 )(y1 − y0 ) − ( x1 − x0 )(y2 − y1 )| p ( x2 − x1 )2 + ( y2 − y1 )2 Hs =
Hsu dd + Hsd du du + dd
(3.8)
(3.9)
where d is the distance from the centroid of a grid cell to a CS; (x0 , y0 ), (x1 , y1 ) and (x2 , y2 ) are the coordinates of the centroid of the grid cell and the two FBPs of the CS, respectively; and the subscripts c, u and d denote the grid cell, upstream and
86 downstream. The streambed properties including the streambed hydraulic conductivity, thickness and top elevation are also specified by users at each SBP. A function of CSR will interpolate these values to each submerged grid cell through the inverse distance method. Because the properties on the bank slope can differ significantly from the main channel of the streambed, using all the SBPs in the property interpolation may remove or weaken the variability among them. Therefore, rather than using all the SBPs, only two nearest upstream SBPs and two nearest downstream SBPs are used in the property interpolation. The stream-aquifer interaction ceases if the stream stage falls to the position equal to or lower than the streambed top elevation. With the flexibility of SBP specification, accordingly, CSR is able to simulate the stream-aquifer interaction where islands or sandbars occur in the middle of the streams or the stream is divided into subchannels in low flow conditions, for example, the Platte River in Nebraska (Smith, 1971; Chen, 2007).
3.3.6
Streamflow routing
The Muskingum-Cunge streamflow routing scheme is applied in CSR to stream segments for predicting the outflow at downstream cross-sections. The MuskingumCunge method is a variant of the Muskingum method developed by Cunge (Cunge, 1969). It is widely used in hydrologic models, such as the SWAT model (Arnold et al., 1998) and the Xinanjing model (Zhao, 1992). Combining the continuity and the storage equations, this method is comparable to diffusion wave routing and provides stable solutions of adequate accuracy by using physically based coefficients and straightforward algorithm (Garbrecht and Brunner, 1991). The Muskingum method is able to account for flood wave downstream movement and attenuation
87 (Fread and Hsu, 1993):
Qd2 = C1 Qu2 + C2 Qu1 + C3 Qd1 + C4 QL
(3.10)
where Qu , Qd and QL denote the upstream flow (inflow), downstream flow (outflow) and lateral flow of the stream segment, respectively; the subscripts 1 and 2 represent the starting and ending time of the routing period; and C1 , C2 , C3 and C4 are the routing coefficients that can be computed by:
(∆t/κ ) − 2χ (∆t/κ ) + 2(1 − χ) (∆t/κ ) + 2χ C2 = (∆t/κ ) + 2(1 − χ) 2(1 − χ) − (∆t/κ ) C3 = (∆t/κ ) + 2(1 − χ) ∆t C4 = (∆t/κ ) + 2(1 − χ)
C1 =
(3.11) (3.12) (3.13) (3.14)
where ∆t is the length of routing period; and κ and χ are the storage coefficient and the weighting coefficient, respectively. To calculate the values of κ and χ, the method developed by Cunge (1969) is used: Ls c
(3.15)
1 q (1 − ) 2 cS0 Ls
(3.16)
κ= χ=
where Ls is the length of the stream segment; c is the flood wave celerity which is approximated by (5(Vu + Vd )/6); Vu and Vd are the upstream and downstream velocity, respectively; q is the unit width flow; and S0 is the average streambed slope of the stream segment. If the calculated value of χ exceeds 0.5, the program
88 will reduce it to 0.5. The value of c is related to the downstream velocity (flow). Additionally, stream-aquifer seepage, dependent on the upstream and downstream stream stages, implicitly affects each other with the downstream flow. Thus Equation 3.10 needs to be solved iteratively by executing the following steps: 1. The initial guess: Qd2 = Qu2 ; 2. Calculate the stream stage Hsd and velocity Vd at DNCS; 3. Calculate the stream-aquifer seepage for each grid cell Qgwc using Equation 3.1 and sum up the total groundwater flow Qgw to or from the current SEG; 4. Reset Qgw = Qin /2 if Qgw < 0 (water flows from the stream to the aquifer) and Qin − Qout < 0, where Qin and Qout are the total inflow and outflow of the SEG, respectively; this step prevents the SEG becoming completely dry; 5. Reset the stream evaporation Qe (this term is always less than zero) to zero if the sum of inflow (positive) and outflow (negative) is less than zero; 6. Calculate the lateral flow QL = Qo + Qe + Qp + Qgw + Qt , where Qo is the overland flow; Qp is the precipitation; and Qt is the tributary flow; 7. Calculate the κ and χ using Equation 3.15 and Equation 3.16; 8. Reset χ to 0.5 if χ > 0.5; 9. Calculate C1 , C2 , C3 and C4 using Equation 3.11 to Equation 3.14; 10. Calculate Qd2 with Equation 3.10; if Qd2 − Q0d2 < TOL (where Q0d2 is the calculated values of Qd2 in the previous iteration and TOL is the flow tolerance,
89 a calculation closure criterion), the calculation is terminated; otherwise start from Step 2 with the new calculated value of Qd2 . When the routing period is relatively long or it is in steady-state conditions, a mass conservation method, which is also used in the SFR packages (Prudic et al., 2004), is applied to simulate flow movement. This method only accounts for mass conservation without flood wave attenuation. Therefore, Step 7 and 8 are skipped, and in Step 9, one obtains:
C1 = 1, C2 = 0, C3 = 0, C4 = 1
3.3.7
(3.17)
Implementation of CSR with MODFLOW
Using a similar program structure as other MODFLOW packages, the CSR package encompasses five main subroutines: AR, RP, AD, FM and BD. The AR subroutine reads the CSR data and allocates spaces for the variables. The RP subroutine reads the data, such as inflows of headwater SEGs and lateral flows of all the SEGs, for each stress period. The AD subroutine resets and initializes the variables at the beginning of a new time step, for example, updating values of the upstream and downstream streamflow. The FM subroutine includes the computation of stream-aquifer seepage and streamflow routing. The exchanged flow is added to the HCOF and RHS terms in Equation 2.14 as follows:
HCOF = HCOF − CDT
Hgw > BOT
(3.18)
RHS = RHS − CDT · Hs
Hgw > BOT
(3.19)
RHS = RHS − CDT( Hs − BOT)
Hgw ≤ BOT
(3.20)
90 where all the terms in the equation are used to describe an individual grid cell that is superimposed by a SEG. The BD subroutine computes the rates and volumes of stream seepages and saves or prints results to specified files.
3.4 Application of the CSR package To test the capability of the CSR package, the proposed approach is applied to simulate the groundwater flow and stream-aquifer interaction of a hypothetical stream-aquifer system (Figure 3.5). In the hypothetical model, the stream network is composed of a diversion channel, a tributary stream channel and a main stream channel. The CS numbers start with 1, 2 and 3 for the main, diversion and tributary channels, respectively. The diversion channel diverts stream water from the main channel to the tributary channel. There are eight SEGs and nine CSs in the main channel, one SEG and two CS in the diversion channel and four SEGs and five CSs in the tributary channel, respectively. The tributary channel joins with the main channel in the middle at SEG [10], which could markedly increase the streamflow at CS [106]. Downstream stages and flows are simulated with the CSR package based on different furthest upstream conditions, which are stream stages or streamflows at CS [101], [201] or [301]. The central part of the streambed of CS [104] is uplifted to form an island in the stream channel when the stream stage is low (Figure 3.5). Table 3.2 shows the properties of each channel. The streambed of the main channel is assumed to have a relatively high hydraulic conductivity to allow sufficient water to flow through for the pumping wells. The single-layer aquifer is divided into 50 rows and 50 columns with a uniform cell size of 100 m in each direction. The thickness of the aquifer is 60 m with a specific yield of 0.2. Hydraulic conductivity is 150 m/d in the vicinity of the main
10 1
91
W1
1
2 10
CS [104]
2
5
3
8
Water Supply Well Stream cross section
5 10
9
10
W5
Stream segment MODFLOW Grid
1
11
W6
12
W7
13
W8
107
GHB
0.5
W4
106
Explanation
0
W3 4 10
4 303
7
302
301
6
202
W2 3 10
304
201
108
2 km
109
Figure 3.5: Schematic of the hypothetical stream-aquifer model. CS [101] to CS [109] represent the main channel; CS [201] to CS [202] represent the diversion channel; and CS [301] to CS [305] represent the tributary channel. stream channel and gradually decreases to 1 m/d towards the model boundaries. There are eight municipal water supply pumping wells evenly distributed along the main channel. Each of the wells keeps a constant pumping rate of 8×104 m3 /d, resulting in a total pumping rate of 6.4×105 m3 /d from the aquifer. The General Head Boundary (GHB) is used at the locations shown in Figure 3.5. Hydraulic heads are set as 52, 50, 48 and 45 m for the GHBs at the northeast side, west side, southwest side and southeast side, respectively. Four steady-state and one transient-state simulations are conducted. The SFR2 package (Niswonger and Prudic, 2005) and the HEC-RAS model (Brunner, 1995) are also selected to compare the simulated results with the CSR package as they are
92 Table 3.2: Stream channel properties in the hypothetical model. Channel Main Length (m) Width (m) Slope (m/m) n Kb (m/d) M (m)
Tributary
5620 2910 250 – 580 117 – 212 0.0015 0.002 0.025 – 0.035 0.025 100 5 3 2
Diversion 1310 20 0.001 0.02 0.01 2
best-known stream routing package in MODFLOW and stream analysis model, respectively. SFR is used in both steady-state and transient-state simulations while HEC-RAS is only used in the transient-state simulation.
3.4.1
Test case: steady-state model
In the steady-state simulations, the inflow of the main channel (at CS [101]) changes between 20 m3 /s and 800 m3 /s respectively as the low flow and high flow conditions . The pumping wells also are deactivated and activated in different cases. By combining different conditions, therefore, there are four scenarios simulated in steady-state conditions: (1) low flow without pumping (LFNP); (2) high flow without (HFNP); (3) low flow with pumping (LFWP) and (4) high flow with pumping (HFWP). The simulated contours of hydraulic heads and the distribution of streamaquifer seepage are shown in Figure 3.6. It is clearly shown that the stream-aquifer seepage varies not only quantitatively but also spatially. In the low flow scenarios, the interaction area of the main channel is reduced to the central part except for CS [104] where the in-stream island is formed and no stream-aquifer interaction occurs. It can also easily be seen from the map that difference of Qgw occurs on a same
93 cross-section from one side to the other, which is more evident in scenarios with pumping. The difference can fundamentally be attributed to the heterogeneity of the streambed geometry and hydraulic properties, which influence the distribution of underlying groundwater hydraulic heads as well.
94 49
W1
W1
49
47
47 W2
W2
45
W3
49
W3
W4
W4
Explanation
Water Supply Well
47
Stream segment Qgw (m3/d)
Water Supply Well
W7
Hgw contour H gw Contour
W6
45
W7
W8
3000
0
(a) Low flow without pumping.
43
Qgw (m3/d)
41 -4000
W6
Stream segment
43
H Hgw contour gw Contour
W5
47
W5
Explanation
-4000
W8
3000
0
(b) High flow without pumping.
43
47
41
4 W1 5
45
47
41
47
W1
41
W2
43
43
W2
3 W3 5
3 W3 9
41
49
37
35
3 W4 3
W4
(c) Low flow with pumping.
35
33
W7 35
Qgw (m3/d)
37
27
35
2200
39
0
29
Qgw (m3/d) -15200
Stream segment Hgw contour H gw Contour
29
Hgw Contour Hgw contour
45
Water Supply Well
31
35
45
Stream segment
31
43
Explanation
W6
39
47
27
31
Water Supply Well
47
33
Explanation
35
33
33
2 W5 9
-15200
0
2200
31 W8
(d) High flow with pumping.
Figure 3.6: Distribution of stream-aquifer seepage and contours of groundwater hydraulic heads near the streams. Hgw and Qgw denote the groundwater hydraulic head and stream-aquifer seepage, respectively. Negative Qgw represents discharge from the stream to the aquifer and vice versa.
95 Table 3.3: Volumetric water budget of steady-state simulations. Positive and negative values refer to gains and losses to the aquifer, respectively. Without Pumping
With Pumping
Low Flow
High Flow
Low Flow
High Flow
CSR Results CSR (m3 /d) -23360.34 3 GHB (m /d) 23376.91 3 Well (m /d) 0.00 Net (m3 /d) 16.57 Discrepancy (%) -0.02
-13344.78 13271.76 0.00 -73.03 -0.07
593768.87 46144.64 -640000.00 -86.49 -0.01
607926.38 32047.15 -640000.00 -26.47 0.00
SFR Results SFR -23551.63 GHB (m3 /d) 23528.86 3 Well (m /d) 0.00 3 Net (m /d) -22.77 Discrepancy (%) -0.01
-14880.88 14613.00 0.00 -267.87 -0.07
592692.04 47193.14 -640000.00 -114.82 -0.02
601681.88 38251.30 -640000.00 -66.83 -0.01
(m3 /d)
Discrepancies of water budget in steady-state simulations are less than or equal to 0.02 % except for the HFNP scenario in which we set a slightly higher tolerance value as the closure criterion of calculation in order to achieve convergence (Table 3.3). However, the discrepancy of the HFNP model is still acceptably low. Table 3.3 also shows that more water is supplied by the stream for pumping in the HFWP scenario than in the LFWP scenario, indicating the model is functioning properly. Besides, the CSR simulated total streambed seepage also agrees with the SFR2 simulated results in all scenarios, because they reproduce matching stream stages and SFR2 uses the actual wet perimeter to calculate the streambed conductance. However, distribution of stream-aquifer seepage is not correctly reproduced with SFR2 due to its limitation in physical representation of wide channels. Figure 3.7 demonstrates the variation of the streamflow and stream stage along the main stream channel. Except for SEG [10], which receives the tributary flow, the streamflow is decreasing in the downstream direction because, the stream is
96 losing water to groundwater pumping wells but not gaining water from other sources (precipitation and overland flow are ignored in the model). 1 2 0 0 5 0 tr ib u ta r y in f lo w
1 0 0 0
4 4 4 0 0
s tr s tr s tr s tr
2 0 0
e a m flo w e a m flo w e a m s ta e a m s ta
( h ig h f lo w ( lo w f lo w s g e ( h ig h f lo g e ( lo w f lo w
1 0 0 0
2 0 0 0
s c e n a r io ) c e n a r io ) w s c e n a r io ) s c e n a r io )
4 2
S tre a m
S tr e a m flo w ( m
6 0 0
s ta g e (m )
4 6
3
/s )
4 8 8 0 0
4 0
0 0
3 0 0 0
D is ta n c e to th e fa r th e s t u p s tr e a m
4 0 0 0
5 0 0 0
6 0 0 0
c r o s s - s e c tio n ( m )
Figure 3.7: Streamflow and stream stage at each cross-section along the main channel (steady-state simulation with groundwater pumping).
3.4.2
Test case: transient-state model
In the transient-state model, the inflow of CS [101] is assumed to change with time as a sine function (Figure 3.8). Due to the short length of the main channel, the time steps therefore are also small to exhibit the flood wave movement and attenuation. In this simulation, the stress periods are uniformly set as 2 minutes and the simulation comprises 144 stress periods. The length of a cycle of the inflow is 72 minutes which includes 36 stress periods. The diverted flow and the inflow at CS[301] are 10 m3 /s and 100 m3 /s, respectively. Streamflows at different CSs are plotted in Figure 3.8. CSs [105] and [106] are the upstream and downstream CSs of SEG [10], respectively. As moving to downstream CSs, the peak flow is attenuated whereas the trough flow is amplified,
97
S tr e a m flo w ( m
3
/s )
1 0 0 0 8 0 0
C S [1 0 1 ]
C S [1 0 5 ]
C S [1 0 6 ]
C S [1 0 9 ]
6 0 0 4 0 0 2 0 0 0 0
5 0
1 0 0
1 5 0
T im e ( m in u te )
2 0 0
2 5 0
Figure 3.8: Hydrographs simulated with the CSR package at different cross-section in the transient-state model. resulting in increasingly small amplitudes. A small increment occurs between CS [105] and [106] due to the tributary inflow for SEG [10]. The outflows of the main channel simulated with different models are also plotted in Figure 3.9. HEC-RAS is used to simulate streamflow in the main channel without the diversion and tributary segments. The simulated flood wave movements (occurrence times of the peak flows) are comparable between CSR and HEC-RAS. The flow simulated with CSR is larger becuase it includes the tributary which is 100 m3 /s. In contrast, SFR2 produces larger lag time for the peak flow with slight attenuation. The peak timing difference between SFR and CSR is approximately one quarter of the inflow cycle period. Because the groundwater seepage is relatively small compared to the streamflow, the difference is induced by the routing methods. HEC-RAS solves Saint Venants equation whereas SFR2 uses a kinematic wave approximation and CSR employs the Muskingum-Cunge scheme which is a finite-difference approximation of the diffusion wave equation. Therefore, CSR and HEC-RAS account for food wave attenuation while SFR2 only accounts for flood wave movement.
98
S tr e a m flo w ( m
3
/s )
1 0 0 0
In f lo w
8 0 0
C S R
S F R 2
H E C -R A S
6 0 0 4 0 0 2 0 0 0 0
5 0
1 0 0
1 5 0
T im e ( m in u te )
2 0 0
2 5 0
Figure 3.9: Comparison of simulated outflow hydrographs at CS [109] using CSR, SFR2 and HEC-RAS.
3.5 Discussions 3.5.1
Advantages
The CSR package uses a cross-section based scheme to represent streams, and couples the Muskingum routing method and a rapid intersection area algorithm with MODFLOW, therefore providing some unique advantages in stream-aquifer simulation. These advantages can be concluded as follows: (1) By using a crosssection based scheme, the stream geometries are physically represented, which overcomes the limitation of other stream packages in representing channels with a width spanning over multiple MODFLOW grid cells. This is also the motivation for developing this package. The package is therefore expected to be used in the circumstance when wide channel streams exist or fine grids are used in the groundwater flow model. (2) It can represent variation of the streambed properties both longitudinally and transversely by streambed points, and account for the effects of variability of streambed thickness, permeability, and roughness on stream-aquifer seepages. Using sufficient streambed points, the package can
99 represent the streambed geometry even of extremely high complexity. (3) This methodology can dynamically predict longitudinal and transverse spatial variability of stream-aquifer interactions, as the streamflow and groundwater levels change. (4) Computational efficiency is improved for stream routing because it is cross-section based rather than grid cell based in which way the actual computational units can be reduced considerably. (5) The package improves ease of data input, especially for complicated stream networks with long stream segments. Rather than assigning the parameters for each cell individually, they are interpolated based on the characteristics of the upstream and downstream cross-sections, which can be relatively easily retrieved using the GIS technology. Besides, the package is independent of the grid discretization that is internally accounted for by the algorithms of CSR. The CSR package can be implemented with varied grid discretization, for example, with local grid refinement.
3.5.2
Assumptions and limitations
The water surface is assumed level across the stream. Stream water flow is perpendicular to the cross-section. The stream-aquifer flux is vertically projected on the horizontal plane whereas the horizontal components of exchanged flux on the sloping streambeds are neglected. The neglected component, however, can play a significant role in the stream-aquifer interactions for streams with deep banks or highly sloping stream bottoms. We suggests that the streambed hydraulic conductivity can be assigned with larger values to practically compensate for ignoring horizontal seepage components. The Muskingum method is used to simulate the downstream movement and peak attenuation of the hydrograph. However, it does not account for backwater
100 effects and it is incapable of simulating very rapidly rising hydrographs in flat channels with slope less than 0.0002 (Garbrecht and Brunner, 1991). The Muskingum method is a numerical approximation of the diffusion wave equation. To minimize the numerical errors, it’s suggested that the Courant number Ls /c∆t is equal to or slightly greater than one (Ponce and Changanti, 1994). The stream segment length Ls is actually fixed after defining the CSs; therefore, the computational time step ∆t needs to be adjusted accordingly based on the flow condition. Another routing method is the simple mass conservation method which assumes the flow instantaneously moves downstream without mass and energy losses. The two methods can make significantly different predictions. Therefore, the mass conservation method is suggested for use only in steady-state conditions, as well as for short stream segments or long computational time steps. The streambed elevation and stream stage at a MODFLOW grid cell are interpolated at the cell centroid based on altitudes of the nearest SBPs and stream stages at the UPCS and DNCS, respectively. When the interpolated stream stage is lower than the streambed elevation, CSR will consider that as a dry stream cell and no stream-aquifer interaction is simulated in this cell. However, a portion of the cell can still be submerged, as the case may be, especially when the grid cell size is relatively large. This type of errors is inherited from numerical difference discretization and cannot be avoided. But it can be reduced by using smaller cell sizes. The current version of CSR does not account for the unsaturated flow when the groundwater level drops below the streambed. This may give unrealistic results of stream-aquifer seepages when the underlying groundwater is considerably deep compared to the streambed position or during the transition from connected to disconnected conditions between surface water and groundwater (Brunner et al.,
101
2011). At the current stage, the package is not compatible with other stream packages though they can run in the same model. It is not allowed to route water between different packages. However, this function has been planned to be implemented in the future.
3.5.3
Potential future improvements
According to the aforementioned limitations, we propose some possible improvements on the CSR package as follows: (1) Apply temporal and spatial refinement to achieve model precision. The local grid renement method (Mehl and Hill, 2005) and variable time step discretization method can be incorporated into the model to reduce numerical errors and instability. (2) Incorporate the capability of simulating unsaturated flow for situations where groundwater falls below the streambed. (3) Incorporate Saint Venant’s equation in streamflow routing to account for backwater effects. (4) Account for the horizontal components of the stream-aquifer seepage on the sloping or vertical banks.
3.6 Conclusions It is very challenging to simulate the streamflow and stream-aquifer routing for wide channels with previous stream packages. In this study, a cross-section based stream routing package is developed for MODFLOW. This approach uses a cross-section based routing scheme, thus being able to effectively represent the variability in the streambed properties such as morphology, permeability and roughness. Accordingly, heterogeneity of the stream-aquifer interactions can be simulated along and across the stream channels. The model was successfully
102 applied to simulate a hypothetical stream-aquifer system in both steady and transient states, as well as with and without groundwater pumping. CSR represents a certain improvement over previous MODFLOW streamflow packages by providing the efficient cross-section based scheme and the unique capability of simulating streambed heterogeneity in longitudinal and transverse directions. This package is applicable to the aquifer systems with wide stream channels, or undergoing stream-aquifer exchange of high heterogeneity. The package is available upon request by contacting the author.
Acknowledgements The research was supported by Lower Platte North Natural Resources District of Nebraska, US National Science Foundation award ID0535255 and China Natural Science Foundation (project no. 41072183 and 41101015).
103
Chapter 4 AN INTEGRATED SURFACE- AND GROUND- WATER MODEL 4.1 Introduction Integrated Water Resources Management (IWRM) approach has been accepted internationally as the way forward for efficient, equitable and sustainable development and management of the world’s limited water resources and for coping with conflicting demands. An integrated hydrologic model is a useful tool for IWRM by projecting the effects of hydrologic processes and human activities on surface water and groundwater interactions, which can substantially influence both the quality and quality of water resources. Over the last two decades, considerable progress in computer and hydrologic modeling techniques has led to the development of physically based, spatially distributed codes that are capable of simulating integrated surface-subsurface hydrological processes. Among them, there are some sophisticated models solving the full three-dimensional Richards equation such as HydroGeoSphere (Therrien et al., 2010), InHM (VanderKwaak, 1999), ParFLOW (Kol-
104 let and Maxwell, 2006), and MODHMS (Panday et al., 2009). However, these models are computationally intensive and highly data-demanding for modeling watershed scale hydrologic dynamics. It involves excessive work in data preparation, model design, calibration and verification, especially for larger scale cases (Partington et al., 2013). In this study, therefore, the Integrated Surface- and Ground- Water Model (ISGWM) is developed by integrating SWAT and MODFLOW, which have been extensively tested and gained wide acceptance in their fields.
4.1.1
SWAT
SWAT is a river basin or watershed scale model developed by the USDA Agricultural Research Service (Arnold et al., 1998; Neitsch et al., 2005). SWAT is capable of predicting the impact of land management practices on water, sediment and agricultural chemical yields in large complex watersheds with varying soils, land uses and management conditions over long periods of time. SWAT is a semidistributed and partially physically-based model. The inputs for SWAT can be easily prepared and its computation is very efficient. The model runs based on daily time step, thus it’s limited to simulate detailed single flood event. SWAT has gained international acceptance in hydrological modeling since its development. It has been proven to be an effective tool for assessing water resource and non-point source pollution problems for a wide range of scales and environmental conditions across the globe. An extensive SWAT literature database is maintained by the Center for Agricultural and Rural Development, Iowa State University (https://www.card.iastate.edu/swat_articles/index.aspx). SWAT has been used to quantify non-point pollution transportation in stream networks (Cheng et al., 2007; Hao et al., 2004; Kalin and Hantush, 2006; Santhi et al., 2001).
105 Mapfumo et al. (2004) applied the model to simulate soil water content in a small watershed under three grazing intensities in Alberta, Canada. Soil water were overestimated in dry soil conditions and underestimated in wet soil conditions. However, the model achieved adequate accuracy in modeling soil water patterns for all three watersheds based on a daily time step. Abbaspour et al. (2009) calibrated a SWAT model in Iran for the period from 1980 to 2002 using daily river discharges and annual wheat yield data at a subbasin level, and assessed the climate impact under the future scenarios downscaled from the predictions of the Canadian Global Coupled Model. They predicted that wet regions of Iran will become wetter whereas dry regions will become drier. However, there are some limitations in SWAT. For example, compared to other components in SWAT, its groundwater flow component is known to be oversimplified. SWAT simulates two aquifers in each hydrologic response unit (HRU), including a shallow unconfined aquifer and a deep confined aquifer (Arnold et al., 1998). The shallow aquifer contributes baseflow to streams within the subbasin in which it is located. However, the groundwater flows and heterogeneities of the aquifer properties, such as conductivity and storage coefficient in the aquifers are ignored by SWAT. The soil water module in SWAT is a simple cascade storage routing model by incorporating a number of parameters. These model parameters, which are difficult to be physically measured, need to be calibrated against observation data, typically streamflow. Even though adjusting a number of lumped parameters can archive high model performances, the individual physical process or the system under study may be inadequately described or even misinterpreted (Srinivasan et al., 2005). Besides, combination of empirical components is prone to the equifinality (non-uniqueness in parameterization) problem, without understanding the true physical meaning of these parameters (Gupta and Sorooshian, 1983;
106 Beven and Freer, 2001; Yang et al., 2008). In the SWAT soil module, it assumes that the soil water is driven only by the gravitational force and water flow occurs only when a layer exceeds its field capacity. Physically, the soil water module does not simulate water movement beyond the soil zone which is usually several-meter deep (Chen et al., 2011). Simulation of deep soil moisture, regardless of the SWAT parameterization chosen for root-water extraction, is limited. Deep percolation, which is considered to become groundwater recharge eventually, is applied to the groundwater table instantly. Furthermore, it’s difficult to implement precise boundary conditions, such as infiltration, surface ponding and groundwater table (Yang et al., 2009). Concurrent hydrologic processes, such as percolation, root uptake and lateral flow are calculated separately. Moreover, another limitation in SWAT is that it cannot account for saturation excess runoff (Easton et al., 2008). Therefore, a more physically-based representation of soil water is needed to simulate the soil hydrologic characteristics and processes.
4.1.2
Integration of SWAT and MODFLOW
There have been a few studies attempting to couple SWAT and MODFLOW. For example, due to the demand of integrated water management and the need of improving groundwater modeling accuracy, an integrated model, referred as SWATMODFLOW (or SWATMOD), has been developed (Perkins and Sophocleous, 1999; Sophocleous et al., 1997). SWATMOD is an interface between SWAT and MODFLOW, which is able to convey data of the water exchange between surface water and groundwater on each HRU. The model was used to simulate the flow of surface water, groundwater, and stream-aquifer interactions. Perkins and Sophocleous (1999)
107 applied the model to analyze the drought impact and examine the relative contributions of stream yield components. Sophocleous et al. (1997) used a two-way coupling conceptualization to modify the interface and successfully applied the model to investigate the irrigation effects on the streamflow and groundwater storage. Galbiati et al. (2006) coupled SWAT, MODFLOW, MT3DMS and the in-stream water quality model QUAL2E in different temporal scales. They successfully evaluated and predicted the surface and subsurface water quality and quantity affected by anthropogenic activities in the Bonello watershed, Italy. Kim et al. (2008) developed the HRU-cell conversion interface, which is capable of effectively simulating the distributed groundwater recharge rate and the groundwater evapotranspiration. The model was successfully tested in the Musimcheon Basin, Korea. The HRU-cell conversion interface requires uniform square sizes for the MODFLOW grids and presume the MODFLOW sizes are smaller than HRUs. In these studies, however, SWAT and MODFLOW were loosely coupled without representing their physical interface, the groundwater table. Therefore, groundwater recharge estimation largely depends on the performance of the runoff model in these study cases. Soil is the inter-medium between SWAT and MODFLOW. Soil moisture is an essential quantity to be represented in dealing with the exchange of energy and water within the soil-vegetation-atmosphere continuum (Chen et al., 2011). Soil water content determines the water supply at the interface between the soil surface and the atmosphere. The groundwater recharge and the generation of surface runoff leading to flood formation are also linked with the soil moisture content. These processes require representations of soil moisture content as a function of time at different scales. Soil moisture is of high spatiotemporal variability due to the combined effects induced by factors such as precipitation, soil texture, drainage, irrigation, flooding and shallow groundwater (Overgaard et al., 2006). Heterogeneity
108 of these individual factors can also increase the spatial variation of the patterns of soil geobiochemical and hydrological processes, controlling a wide array of ecological, hydrological, geotechnical, and meteorological processes. The role of these factors in controlling the spatial patterns and temporal dynamics is often not well known (Korres et al., 2015). Chen and Hu (2004) proved that considering groundwater effects in the soil water model increased evapotranspiration by over 20%, and further pointed out that variations in the groundwater table can create an additional spatial-temporal variability of soil moisture and surface water flux. Thus, it is of critical importance to develop a soil water module that can physically represent the soil-vegetation-atmosphere continuum and links SWAT and MODFLOW in ISGWM.
4.1.3
Objectives
ISGWM is developed by coupling SWAT and MODFLOW with a soil water module, which is originally developed by Ross (2003). The Ross’ solution is selected because of its high computational efficiency, mass conservativeness and most importantly, its accuracy in representing the physical soil water processes. The soil water module can calculate the exchange flows between SWAT and MODFLOW as shown in Figure 4.1. ISGWM is applied to the Johnson Creek watershed located within the Lower Platte basin, Nebraska. The model accounts for the climate condition, land use, as well as watershed and hydrogeologic properties. The land surface processes, groundwater flow, and their interactions are simulated simultaneously. Irrigation pumping is estimated based on the soil water deficit or crop water stress.
109
Exchange flow SWAT
Direct runoff ET
Infiltration
SWM
STREAM Groundwater ET
Recharge
MODFLOW
Stream-aquifer seepage
Figure 4.1: Exchange flows between the models.
4.2 Methods 4.2.1
Mathematical description of the soil water module
4.2.1.1
Richards equation
The Richards equation, incorporating Darcy’s law and the continuity equation, is currently the most accurate and reliable approach to simulate soil water movement in variably saturated conditions. The Richards equation can be expressed in several forms using different dependent variables. These equations are generally solved numerically to model the soil water movement and solute transport. However, its instability (non-convergence), substantial computational time and the mass balance errors produced by the numerical solution have obstructed its use in distributed hydrologic models. To avoid these obstacles in ISGWM, therefore, we implemented
110 a highly efficient soil water module (SWM) originally developed by Ross (2003). In the Ross’ solution, the Kirchhoff potential (φ, [L2 T −1 ]) is introduced as a variable to represent the soil hydraulic condition:
φ=
Z h −∞
K (h)dh
(4.1)
where h is the soil water pressure [L] and K is hydraulic conductivity [LT −1 ]. The Kirchhoff transform is commonly used to linearize the differential equation. Likewise, with Equation 4.1, we can reduce the original Darcian flux equation that has a water pressure head dependent hydraulic conductivity K into a linear form and thus the Richards equation can be written as:
q=−
∂φ +K ∂z
∂q ∂θ = − − Ss ∂t ∂z
(4.2)
(4.3)
where q is the Darcian flux [LT −1 ]; θ is the soil water content [L3 L−3 ]; t is time [T]; z is the vertical coordinate measured positive downwards [L]; Ss is the net flow of source and sink terms [T −1 ]. Ss can be composed by root water uptake, lateral drainage, septic tank leakage or tile drainage. While the focus here is on a one-dimensional representation of soil water, the air phase flow and solute and heat transport are not considered in this study. However, variation of temperature, fate and transport of nutrients such as nitrogen or phosphorus are simulated in SWAT. Nitrogen, for example, can be transported by the upward convection flow in the top layer due to soil evaporation. Soil temperature is also estimated at each soil layer at different depths. While some of the variables can influence the soil water directly, for example blocking induced
111 by frost conditions, soil hydrological dynamics simulated by SWM feeds back with the soil moisture and flux conditions needed for modeling other processes in SWAT. 4.2.1.2
Soil constitutive functions
The constitutive relationships among the soil water content θ, pressure head h, and hydraulic conductivity K can be described by the Brooks-Corey or a modified Van-Genuchten models (Crevoisier et al., 2009). In SWM, the Brooks-Corey model (Brooks and Corey, 1964) is used: θ − θr = S= θ s − θr
h hb
−λ (4.4)
K = Ks Sη
(4.5)
where S is the effective saturation; θr and θs are the residual and saturated water content, respectively; hb is the bubbling or air entry pressure head [L]; Ks is the saturated hydraulic conductivity [LT −1 ]; λ and η are the coefficients describing the shape of the soil retention curve. S is equal to one when h > he . According to Equation 4.1, we can express φ by using S or h as dependents: where φb = Ks hb /(1 − λη ). In frost conditions, the soil permeability drops dramatically due to blocking of immobile frozen soil water (Burt and Williams, 1976). The frost hydraulic conductivity K f [LT −1 ] is adjusted as (Maidment, 1993): K f = K · min(0.001, 2.0 − 1.9
where θ f is the water content at the field capacity.
θ ) θf
(4.6)
112
4.2.2
Numerical implementation of SWM
Because the unsaturated flow equation is highly non-linear, it is very difficult to find the analytical solution except for restrictively simplified cases with a number of ideal assumptions (Celia et al., 1990; Feddes et al., 1988). Thus, the parabolic partial differential equation is typically numerically solved in practical use (Browne ˇ unek et al., 2008; Sim et al., 2009; van Dam and Feddes, 2000). The standard Picard ˚ and Newton-Raphson linearization methods, which are commonly used in other unsaturated flow models, need to perform the solving processes iteratively until convergence is achieved. The iterative calculation often takes an undesirably long CPU time, especially in watershed models where it needs to be applied to over a large number of soil columns. Therefore, SWM is integrated in SWAT because it is a fast solution based on a non-iterative, adaptive time stepping and mass conservative numerical solution of the Richards equation. Compared with other models, it reduces the computational time significantly and meanwhile accounts for complicated boundary conditions (Varado et al., 2006; Crevoisier et al., 2009). 4.2.2.1
Soil profile
A soil profile is spatially discretized by N layers numbered downwards from the top of the soil profile. The mass balance equation of a soil layer i can be expressed as: ∆zi ∆Si (θsi − θri ) = qiσ−1 − qiσ − (Sri + Sli )∆zi ∆t ! 0 0 ∂q ∂q i i ∆Si + ∆Si+1 qiσ = q0i + σ ∂Si ∂Si+1 q0i
h i φi0 − φi0+1 0 0 = − ωKi + (1 − ω )Ki+1 (∆zi + ∆zi+1 )/2
(4.7) (4.8)
(4.9)
113 where the superscript 0 denotes the beginning of the time step; the superscript σ denotes a fraction the time step; the subscript i denotes the layer number; ∆z is the layer thickness [L]; ∆t is the length of the time step [T]; ω is the weighting factor to calculate the interface hydraulic conductivity; Sr and Sl are the root uptake and lateral outflow, respectively [T −1 ]. The solution of Equation 4.7 is perfect when the exact values of σ and q can be computed. However, it is unachievable or too difficult. Therefore, the value of σ presumed to be 1 if a saturated layer occurs or 0.5 for other cases. Equation 4.8 is actually a first order Taylor polynomial using the variables at the beginning of the time step. When Equation 4.8 is applied to the saturated layer, S needs to be replaced by φ for that layer. The value of ω is updated at each interface and time step by satisfying that the flow calculated with Equation 4.9 is zero when the two layers are hydrostatic. Combining Equation 4.7 and Equation 4.8 of each layer yields a tridiagonal set of linear equations with unknowns being ∆Si or ∆φi : ai ∆Si−1 + bi ∆Si + ci ∆Si+1 = di
(4.10)
where ai , bi , ci and di are the coefficients and can be calculated as follows: ∂qi−1 0 ai = − ∂Si−1 ∆zi (θsi − θri ) ∂qi−1 0 ∂qi 0 bi = − + σ∆t ∂Si ∂Si ∂qi 0 ci = ∂S
(4.11) (4.12) (4.13)
i +1
di =
q0i−1 − q0i − (Sri + Sli )∆zi σ
(4.14)
114 4.2.2.2
Source and sink
Soil water uptake (transpiration) by plants is one of the most crucial processes influencing the hydrological and ecological dynamics in soil. On the other hand, availability of water and nutrient largely constrained plant growth. A simplified version of the Environmental Policy Integrated Climate (EPIC) model is embedded in SWAT to simulate plant growth under the effects of climate, water, and nutrient. Therefore, soil moisture simulated with SWM can in turn improve the plant growth model by providing accurate estimation of available water for plant growth use. Root water uptake is strongly correlated with distribution of the surface area and mass of plant roots (Feddes et al., 2001). The actual water uptake by root systems can be simulated using potential transpiration, water availability factor, root mass distribution and maximum root depth. In ISGWM, the cumulative potential water uptake rate wr,cum [L/T] is estimated as:
wr,cum
Et z = 1 − exp(− β ) [1 − exp(− β)] zr
(4.15)
where Et is potential transpiration rate [L/T]; β is the root distribution factor; and zr is the root depth [L] that is determined by SWAT based on the type of the vegetation and its growing stage. The plant uptake condition will change from the demand-driven state to the supply-dependent state (Raats, 2007) when root water uptake is limited by soil water. In the supply-dependent state, water stress hinders plant roots in creating hydraulic gradients at the soil-root interface. Therefore, the actual amount of water
115 extracted by plant roots also depends on soil water content:
Sr =
w
(θ ≥ 0.25θ a )
r
wr exp 5
θ −1 0.25θ a
(4.16)
(θ < 0.25θ a )
where wr is root uptake rate wr at a soil node [1/T]; and θ a is the available water content. θ a = θ f − θw , where θ f and θw are the soil water contents at the field capacity and permanent wilting point, respectively. In most cases of using the one dimensional Richards equation model, the lateral flow is neglected because the soil column is considered as a closed system horizontally. In watershed hydrology, subsurface flow may contribute a large portion of the total runoff. A kinematic storage model is employed in SWAT to account for lateral flow in the unsaturated zone (Sloan and Moore, 1984). Lateral flows occur only when the volumetric water content is greater than the field capacity: Sl =
2K (θ − θ f )S0 (θs − θ f ) Ls
(4.17)
where S0 is the slope of the soil layer which is assumed as the land surface slope, [LL−1 ]; and Ls is the hill slope length [L]. 4.2.2.3
Top boundary condition
The flux across the surface-subsurface interface can be upward (soil evaporation or exfiltration) or downward (infiltration). Depending on land surface and atmospheric conditions, as well as the ability of soil to transmit water, the exchange processes occurring at the top boundary are extremely complex (Feddes et al., 2001). The hydrological processes simulated by ISGWM occurring at the soil surface
116 include: overland runoff, infiltration and soil evaporation. The potential net flux qnet [L/T] asserted at the soil surface can be described as:
qnet = qon + I + R + Sm − Es
(4.18)
where qon is the flux of run-on from the upstream subbasins; I, R and Sm are the rates of irrigation, net rainfall and snow melt reaching the soil surface, respectively [L/T]; and Es is the potential soil evaporation rate [L/T]. Partitioning runoff and infiltration from rainfall is the core task of hydrologic modeling. Besides the vegetation and atmospheric conditions, runoff/infiltration processes are profoundly influenced by temporal dynamics of near-surface soil moisture conditions. While runoff/infiltration is simulated through the SCS Curve Number or Green-Ampt methods in SWAT, SWM simulates these processes through specifying appropriate boundary conditions. To represent the runoff generation process, SWM formulates the top boundary condition under the ponding or non-ponding conditions. The soil water pressure heads in the top soil are highly varying. In temperate weather and moisture conditions, the potential net flux will not exceed the water transmitting ability of the top soil. In such cases, the flux boundary condition is applied to the top of the soil profile by substituting the upper flow term in Equation 4.14: q00 = qnet
(4.19)
When infiltration predominates and exceeds the infiltration capacity of the top soil, it will prompt water ponding on the soil surface. The depth of ponding water plays an important role in estimation of both the infiltration rate and overland
117 flow (Wallach et al., 1997). Warrick et al. (2005) concluded that the ponding depth has a significant effect on infiltration through modeling infiltration under variable ponding conditions. Philip (1993) evaluated the quantitative effects of ponding depth on infiltration and concluded that ponding can increase the amount of infiltration. Water ponding can be triggered by intense rainfall or during the late stage of irrigation. In such cases, a virtual node is created at the soil surface and it is assumed that the pressure head of the node is equal to the depth of ponding water. Using the falling head method, the variable-ponding-depth boundary condition can be explicitly specified on the basis of mass conservation: ∆h0 = qnet − qσs − q0σ ∆t
(4.20)
where qs is the surface runoff [L/T]; and q0 is the infiltration rate [L/T]. When the Green-Ampt method is used in SWAT, ponding water is entirely converted to overland runoff at the end of each time step after the infiltration calculation. As a result, the partition of rainfall can be affected by the length of the time step. Longer time steps will correspond to longer ponding time, and thereafter increase infiltration. In SWM, the amount of ponding water converted to runoff is calculated with the Manning equation. As shown in Figure 4.2, the overland flow is considered as shallow flow through a very wide rectangular channel and thus hydraulic radius can be approximated as the ponding depth. The flow velocity vov [L] can be estimated as:
√ vov =
S0 ( h 0 ) κ −1 ns
(4.21)
where ns is a resistance factor that depends on the cover of the land surface; and κ
118 is a constant depending on the flow conditions (Table 4.1). While overland flow typically behaves as turbulent flow, κ is assumed to be 5/3 in SWM. Assuming water moves uniformly in the downslope direction, the water removed from (shaded area in Figure 4.2) the slope plan is:
√ vov h0 S0 qs = = ( h0 )κ Ls ns Ls
(4.22)
The surface runoff will be added to the SFR package as overland runoff in MODFLOW and the rest of water will keep ponding on the soil surface in the next time step. Under the ponding condition, the flow equation for the top node is written
S0
Figure 4.2: Uniform movement of overland flow. as follows:
1 dqs ∂q + + 0 σ∆t dh0 ∂h0
∆h0 +
qnet − q0s − q00 ∂q0 ∆S1 = ∂S1 σ
(4.23)
Equation 4.23 is added to the linear equation series and thereby runoff calculation is included in the soil water model. To account for the vegetation and slope effects on the runoff processes, a skin layer is added to the land surface. The land surface effects are justified by the hydraulic conductivity Kq [L/T]) and thickness Mq [L] of the skin layer. Mq is also the minimum ponding depth with which runoff can occur.
119 Table 4.1: Values of κ used in overland flow estimation Turbulence 100% κ
75% 50% 25%
0%
5/3 6/3 7/3 8/3 Source: (Ponce, 1989)
9/3
Based on Equation 4.22, Table 4.2 is an example showing the calculated time required to convert all ponding water to overland runoff on a hypothetical HRU. If the ponding depth is larger than 10 mm, ponding water will be converted to runoff within one hour. However, as the ponding depth decreases, the time increases for ponding water “moving out” from the HRU. Table 4.2: Overland runoff qs and time ts to entirely convert to runoff for different ponding depth h0 h0 (mm)
qs (mm/hr)
ts (hr)
0.1 0.007 13.68 1 0.339 2.95 10 15.754 0.63 100 731.239 0.14 Ls = 100 m; S0 = 0.15; ns = 0.15 SWAT partitions total evapotranspiration into canopy evaporation, plant transpiration, sublimation and soil evaporation. Soil evaporation is often divided into two distinct stages. At the first stage of soil evaporation, soil is relatively wet and thereby soil evaporation is determined by the surface energy balance and mass-transfer conditions such as wind speed and humidity, which is so-called the atmosphere-controlled evaporation (or energy limited) stage. SWAT estimates the potential evaporation rate of soil with accounting for the atmospheric conditions on a daily basis. When soil moisture depletes and the pressure head of the top soil layer drops below a threshold value, it suggests that soil evaporation is constrained by soil moisture, which is so-called the soil-controlled (or soil limited) stage.
120 During the soil-controlled stage, suggested by Feddes et al. (1988), the threshold pressure head value hem [L] can be specified as:
hem = h a =
RT ln( Hr ) Mg
(4.24)
where M is the molecular weight of water [Mmol −1 ]; g is the gravitational acceleration [LT −2 ]; R is the gas constant [Jmol −1 K −1 ] and Hr is the relative humidity. Allen et al. (2005) suggested that the soil evaporation limiting threshold can be adjusted according to soil and climate condition and soil water is only evaporable when θ > 0.5θw . SWM provides two options to define this threshold at the top of the soil profile: (1) automatically computed by the code using Equation 4.24 based on the atmospheric condition; or (2) specified with a soil water pressure. Before each time step, the top boundary type is determined based on h00 and qnet (Figure 4.3). If ponding water occurs at the beginning of the time step, the falling head boundary is applied to the soil surface. Otherwise, a maximum flux q00 is estimated by assuming the the pressure head of the top node is equal to the minimum depth for overland runoff Mq or the evaporation threshold hem , depending on the direction of qnet . If the absolute value of q00 is smaller than qnet , it indicates that the potential net flux qnet exceeds the soil transmitting capacity. In this case, the flux q0 is set as q00 . Otherwise, a flux boundary is used with q0 being equal to qnet . 4.2.2.4
Bottom boundary condition
The bottom boundary is relatively simple by using the groundwater table as the specific head boundary. The flow across the soil bottom is the groundwater net recharge. Because the groundwater head at the beginning of a MODFLOW time
121 Yes
No
h0 > 0 Yes
h0’ = Mq
h0’ = hem
q0’ < qnet
q0’ < qnet
Yes Ponding
No
qnet > 0
No
Initial Ponding
Yes
Flux Boundary
No Restricted ET
Figure 4.3: Top boundary condition types used for SWM. step is used as the boundary, the SWM only needs one run for each MODFLOW time step. 4.2.2.5
Time discretization in SWM
In general, larger time steps can reduce the computer running time but also produce larger numerical errors. In SWM, an automatic time-stepping scheme is used to balance the numerical errors and the computational time. The time-step size is adjusted to maintain two criteria: 1. For all the layers, ∆Si needs to be smaller than a maximum change of the effective saturation ∆Smax specified by the modelers. Therefore, the initial time step length ∆t0 [T] can be evaluated as such: ∆Smax ∆t0 = q i −1 − q i 0 ∆z (θ − θ ) i si ri max
(4.25)
2. The saturation state of each layer cannot change during a time step because the choice of equations is highly dependent on the saturation state. When a
122 saturation sate change is identified at the end of a time step, the time step will be proportionally reduced.
4.2.3
Surface runoff
The generation of surface runoff is simulated by SWM. After that, SWAT calculates the flow available to streamflow routing after accounting for surface runoff lags due to a time of concentration greater than one day in large subbasins. The detail of determination of runoff lag is well documented by Neitsch et al. (2005). Streams are simulated by the Streamflow-Routing Package in MODFLOW (Prudic et al., 2004). Stream water is routed using a kinematic-wave method. The overall water budget of a stream segment can be expressed as:
Qsri + Qtrb + Qro + Q ppt − Qbed = Qsro + Qdiv + Qet
(4.26)
where Qsri is upstream flow or flows from such as lake or model inlets [L3 /T]; Qtrb is the tributary flow [L3 /T]; Qro is the overland runoff [L3 /T]; Q ppt is precipitation that falls directly on a reach [L3 /T]; Qbed is the streambed leakage; Qsro is downstream flow out from the segment [L3 /T]; Qdiv is the specified diversions from the segment [L3 /T]; and Qet evapotranspiration from the segment [L3 /T]. Surface runoff, which is calculated by SWAT for each subbasin, is uniformly distributed to the stream segments as Qro based on the stream lengths. One subbasin must have at least one stream segment. However, it is allowed to designate the subbasin to more than one segments when there are multiple stream segments located within one subbasin.
123
4.2.4
Irrigation requirement
Crop irrigation requirement is estimated based on the soil water deficit or plant water stress thresholds. Soil water deficit is the difference of the required and actual soil moisture. Likewise, plant water stress is the difference of the potential and actual transpiration. At the beginning of the soil water model time step, if there is an irrigation well located in the HRU and the crop irrigation requirement is not met, the model will trigger the pumping well. The total pumping rate of an individual well is:
Qw = min[( Ic /E f − R − h0 /∆t) A I , Cw ]
(4.27)
where Qw is the actual pumping rate of a irrigation well, [L3 /T]; Ic is the net irrigation requirement estimated based on the soil water deficit or plant water stress thresholds or the maximum value of these two, [L/T]; A I is the irrigated area, [L2 ]; E f is the irrigation efficiency; Cw is the pumping capacity of the well. The min function assures that the estimated actual pumping rate is not lager than the pumping capacity. Equation 4.27 accounts for the rainfall and ponding water. Besides the well pumping capacity and irrigation efficiency, the model also determines the status of the well during the simulation based on the well start and abandon times. By doing this, the model can consider the conversion from dryland to irrigated land. The variables needed for an irrigation well are explained in Appendix C. Irrigation deficit is estimated for the HRU designated to the well. However, the locations of HRU and the MODFLOW cell can be different. In such cases, water is transferred from one HRU to another. At the beginning of each SWM step, Qw is estimated based on the crop or soil conditions. IRRDT is used to
124 convert the soil water deficit into flow rate. The smaller the IRRDT is, the more intense is the irrigation rate. Irrigation water is pumped from the groundwater model cell designated with the well. Soil evaporation deficit will be first fulfilled because irrigation water is uniformly applied to the surface of the HRU. Vegetation interception is ignored in the irrigation simulation.
4.2.5
Model interfacing
4.2.5.1
Temporal and spatial discretization
The aquifer is discretized into rectangular blocks in MODFLOW while the watershed is discretized into HRUs in SWAT. Time discretization in MODFLOW is arbitrary though it usually runs at larger time steps such as weekly or monthly for long term water resource assessment purposes. SWAT, on the other hand, operates on a daily time step. SWM employs an auto time stepping scheme whose time step is even smaller. Therefore, it is very important to develop a model interface that correctly conveys the values between models to account for the difference of their spatial and temporal discretization. Overlaying the MODFLOW grid and SWAT HRUs can generate a layer of interface fractions (Figure 4.4). Each fraction has a MODFLOW cell index and a HRU index. Hence, specification of each of the interface fractions connects SWAT and MODFLOW. Scales of SWAT HRUs and MODFLOW cells can be arbitrary. The flow and groundwater head are calculated using the weighted average method. The weights are automatically determined by ISGWM according to the interface fraction areas. For larger scale models, subbasins are commonly represented with fewer HRUs or even one single HRU to maintain efficiency. In such cases, those MODFLOW cells, not covered by the effective HRUs, are conected to the subbasin
125 instead. The average recharge of the subbasin will be used for those grid cells.
Figure 4.4: Interface fractions generated by overlaying MODFLOW grid and SWAT HRUs. In ISGWM, SWAT and SWM become a package for MODFLOW, which is referred to as Unsaturated-zone and Surface Water package (USW). At the beginning of a MODFLOW time step, the depth to groundwater table of each HRU is updated. Before each MODFLOW time step, ISGWM determines whether a USW run is needed. If the MODFLOW time step is larger than one day, the USW will be executed for the number of days in the MODFLOW time step. If the MODFLOW time step does not start at the beginning of a day, the flows from the last day of the previous MODFLOW time step, such as recharge and surface runoff, will be used for the rest of the day and accumulated for the MODFLOW time step. The net recharge and well pumping calculated by USW are used as additional flows in
126 the MODFLOW recharge and well packages, respectively. Thus, ISGWM requires activation of these two packages in MODFLOW during the simulation. 4.2.5.2
Initial conditions
Unlike groundwater flow model, the initial conditions of a watershed model is very difficult to define because of the challenge of measuring the large number of variables, such as stream stages, soil water contents, snow storages and vegetation conditions. Therefore, watershed models commonly need a warm up period to self-accommodate the hydrologic conditions. In ISGWM, the warm-up period is defined by specifying ”NYSKIP” in the SWAT basin file. During the warm up period, the SWAT groundwater module (GWMOD) and streamflow routing module, which will be replaced by MODFLOW, are activated. The groundwater table of the shallow aquifer in GWMOD serves as the lower boundary condition for SWM. On the other hand, MODFLOW is disabled during the warm-up period. 4.2.5.3
Input and output files
During an ISGWM simulation run, the net recharge and irrigation pumping flow are written to disk following the MODFLOW file format. Therefore, when the surface water model is calibrated, these files can be used for groundwater simulation in MODFLOW without running SWAT and SWM. One of the most advantageous features of SWAT and MODFLOW is that there are various tools to pre- and post- process the input and output files of these two models. Therefore, ISGWM directly reads the model files which are set up individually for SWAT and MODFLOW. The extra information needed by ISGWM are the interface fractions, soil properties for SWM and irrigation well information. The interface fractions can be created through the ”Split” function in GIS analysis
127 tools. Soil properties are estimated based on their textures which are specified in the SWAT soil file. A R scipt converting the SWAT soil file to SWM file is provided in Appendix D. The data input instructions for ISGWM are also provided in Appendix C.
4.3 Application of ISGWM 4.3.1
Study Area
To demonstrate the applicability of ISGWM, it is applied to the Johnson Creek watershed (JCW) located within the Lower Platte basin, Nebraska. The Johnson Creek is a tributary of the Silver Creek which joins with the main channel of the Platte River. JCW drains an area of about 56 km2 into the Platte River (Figure 4.5). The land surface elevation decreases southwestward from 414 to 327 m above the sea level with an average slope of 1.3%. The daily stream discharge is measured by the USGS gage station 06804900 located at the watershed outlet. Land use data (Figure 4.6) is obtained from the Center for Advanced Land Management Information Technologies (CALMIT) at University of NebraskaLincoln. As shown in Table 4.3, the primary land use type in the Lower Platte River basin is agricultural farmlands (80.6%). Among the crop types, nearly 42% are corns. Irrigated lands are mainly distributed in the low land which is on the west side of the stream. There are 70 irrigation wells located in JCW for the irrigated farms covering over 40% of the watershed area. Irrigation well properties are obtained from the Nebraska State Registered well database. The Soil Survey Geographic (SSURGO) Database is used to retrieve soil information for JCW. There are over 1000 soil map units in this area according to the
128
A A
¡À
A
A A
A A A
Legend City
d n
A
A
YutanWell Irrigation
Subbasin
A
AA
A A
Stream Gauge Stream
A
A
Weather Station
A
A
MeadA
Elevation
A
Value
A
A
A
A A
A A
A
A A
A
A
A
A
A AA
Low : 327.659
A
A
High : 414.677
A
A
A A
A
A
A
A
A
A
A
Nebraska
A
A
A
A
A A
A A
A
0
1
2
4 Kilometers
A
A
A MEADAGROFARM
d
A
6804900
n
Figure 4.5: Location and topography of the Johnson Creek watershed. Table 4.3: Major land use types in JCW Land use type
Percentage
Irrigated Corn Irrigated Soybeans Range, Pasture, Grass Dryland Corn Dryland Soybeans
23.1% 16.9% 15.4% 18.8% 16.0%
SSURGO map. After aggregating SSURGO Data they can be categorized into four USDA soil classes based on the texture of the whole soil profile (Figure 4.7a). The water storage capacity (θ · ∆z) of each soil layer is also plotted in Figure 4.7b. The parameters for the Brook-Corey model are calculated using the pedo-
129 Legend Land use type Dryland Alfalfa Dryland Corn Dryland Small Grains Dryland Soybeans Irrigated Alfalfa Irrigated Corn Irrigated Small Grains Irrigated Sorghum (Milo, Sudan) Irrigated Soybeans Open Water Other Agricultural Land Range, Pasture, Grass Riparian Forest and Woodlands
0
1
2
Roads
4 Kilometers
Summer Fallow Urban Land
Figure 4.6: Land use of the Johnson Creek watershed.
Legend Soil USDA Texture ClLo SiCl SiClLo Soil Water depth, cm 10
SiLo
4.7
90 20
7.9 14.7
40
90.1 Cl
50
SiClLo
80
SaClLo
µm
ClLo 30
70
SaCl
40
50 2−
SiCl
ilt
50
]S [%
[%
60
60
]C
26.8
70
lay
0−
2µ
m
30
80
20
90
Lo SiLo
0
1
2
4 Kilometers
SaLo
10 Sa
LoSa
Si 10
20
30
40
50
60
70
80
90
[%] Sand 50−2000 µm
(a) Soil distribution.
(b) Soil texture classes and the soil water storage capacity.
Figure 4.7: Soil distribution and textures of JCW.
130 transfer functions (PTF) based on the texture and the porosity (Maidment, 1993): θr = − 0.0182482 + (0.00087269SD ) + (0.00513488CL ) + (0.02939286θs )
− (0.00015395CL 2 ) − (0.0010827SD θs ) − (0.00018233CL 2 θs2 )
(4.28)
+ (0.00030703CL 2 θs ) − (0.0023584CL θs2 ) λ = exp(−0.7842831 + (0.0177544SD ) − (1.062498θs ) − (0.00005304S2D )
− (0.00273493CL2 ) + (1.11134946θs2 ) − (0.03088295SD θs ) + (0.00026587S2D θs2 )
(4.29)
− (0.00610522CL2 θs2 ) − (0.00000235S2D CL ) + (0.00798746CL2 θs ) − (0.00674491CL θs2 )) he = exp(5.3396738 + (0.1845038CL ) − (2.48394546θs ) − (0.00213853CL2 )
− (0.04356349SD θs ) − (0.61745089CL θs ) + (0.00143598S2D θs2 ) − (0.00855375CL2 θs2 ) − (0.00001282S2D CL ) + (0.00895359CL2 θs )
(4.30)
− (0.00072472S2D θs ) + (0.0000054CL2 SD ) + (0.5002806θs2 CL )) η =2/λ + 3
(4.31)
where SD and CL are the relative abundances of sand and clay, respectively. Figure 4.8 shows the comparison of the water contents at -1/3 and -15 bar between the BC model and measurement. The parameters used for the BC model is estimated using PTF. The determination coefficient is over 0.85, indicating the BC model with the parameters estimated with PTF can be used to describe the relationship of the soil properties in this area. As shown in Figure 4.5, there is one weather station, which belongs to the Automated Weather Data Network (AWDN), operated by the High Plains Reginal Climate Center since 1994. The daily solar radiation, wind speed, precipitation, temperature and relative humidity are measured at the station. These are also the climate variables needed in SWAT for estimating potential evapotranspiration using the Penman-Monteith method. The mean annual precipitation of JCW
131
●●
Pressure 0.3
● ●
1/3 Bar
● ●
● ●●
15 Bar
●●
● ●●● ●● ● ● ● ●● ●
●
● ● ● ●
●
●
Simulated
● ● ● ● ●● ● ●● ● ●●● ●
0.2
● ● ● ● ●● ● ● ●● ● ● ● ●● ● ● ●●
0.1
●●
R2 = 0.853 0.0 0.0
0.1
0.2
0.3
Measured
Figure 4.8: Calculated and measured θ of soil in JCW. between 1994 and 2009 is 685 mm.
4.3.2
Model Description
The input files for SWAT is generated using the tool ArcSWAT (Arnold et al., 1998). JWC is subdivided into nine subwatersheds and 251 HRUs. In this exercise, all the HRUs generated by ArcSWAT are included in the simulation. The upper two meter of the soil profile is uniformly discretized into 0.25 m layers while the lower is 0.5 m. The hydrogeoloic properties, including horizontal hydraulic conductivity (Kh ) and specific yield (Sy ), are acquired from the Lower Platte River basin model (Chapter 2). In this exercise, the four-layer aquifer system is reduced to one layer to avoid dry cell issues in MODFLOW. The new parameter values are recalculated using the harmonic mean weighted by layer thicknesses. The vertical hydraulic
132 conductivity (Kv ) is used as the saturated hydraulic conductivity for the lower layers of the soil profile. Meanwhile, Sy is used to estimate θs = θr + Sy , such that continuity between the unsaturated zone and unconfined aquifers is maintained. The MODFLOW grid is regenerated using a uniform grid with a cell dimension of 170 m by 170 m, which is consistent with the resolution of the DEM used for watershed delineation with ArcSWAT. Therefore, each MODFLOW cell coincides with the HRU boundary (Figure 4.9). Initial heads are generated by a steady-state run by specifying constant heads at some groundwater observation well locations. In the transient run, general head boundaries are used at the north and south model boundaries. Fourteen stream segments compose the stream network, with two segments in each headwater subwatershed and one segment in other individual subwatershed. In the headwater subwatershed, surface runoff only contributes to the upstream segment. Streambed elevations are determined using the lowest values of a 30meter DEM within each stream cell and are verified to maintain the down gradient in the downstream direction. The simulated period covers from 1994 to 2009 with the first four years being used as the warm up period for SWAT. In MODFLOW, the time spans from 1998 to 2009 is discretized into 146 monthly stress periods with stress period length assigned following the calendar months. Each stress period is further equally subdivided into two time steps. Therefore, in some months, the second MODFLOW time step starts at the 15.5th day. In such cases, the flows simulated with SWAT in day 16 are used for MODFLOW for the rest 0.5 day. And SWAT runs from the 17th day.
133
Legend Stream
2
1
HRU
MODFLOW grid cell Initial head
6
7
0.000000 0.000001 - 334.954000 334.954001 - 340.852000 340.852001 - 347.123000
8
9
347.123001 - 351.456000 351.456001 - 354.382000
10
354.382001 - 357.995000 357.995001 - 362.636000
11
3
362.636001 - 366.668000
12
1
2
4 Kilometers
14
0
Figure 4.9: Discretization of SWAT HRUs and MODFLOW grid cells for JCW. The numbers along the streams are the segment numbers.
4.4 Simulation results of ISGWM The purpose of this modeling exercise is to demonstrate the applicability of ISGWM. Therefore, it is calibrated only with the monthly streamflow data at the stream gauge station (Figure 4.5). The soil hydraulic conductivity and porosity are manually calibrated to match the simulated streamflow with the measured values. The simulated and measured streamflows are plotted in Figure 4.10. The seasonal variation in streamflow is represented in the model. Overall baseflow contributes over 90% of the total streamflow. The streamflow is slightly underpredicted with a R2 of 0.6. The underestimate of the streamflow may be caused by the streambed hydraulic conductivity which is a difficult term to measure. Net recharge is the net flow across the interface between the unsaturated zone
134 and saturated zone. It is the most important variable simulated in ISGWM. As shown in Figure 4.11, the ratio of the annual recharge and precipitation ranges from 2.42% to 23.36%. Furthermore, the ratio is larger in wet years than dry years. This can be attributed to the better permeability of soil in the wet condition.
06804900
Discharge (m3s)
1.00
Model Observed Simulated
0.01
1998
2000
2002
2004
2006
2008
2010
Date
Figure 4.10: Simulated versus observed monthly streamflows. The distribution of mean monthly net recharge is plotted in Figure 4.12. It is obvious that recharge rate is higher between May and July. In August, net recharge becomes negative in the irrigated area, which may be caused by the net water requirement for ET reaching maximum. Soil water is depleted because of the increasing water use requirement for crops. Irrigation well pumpages are estimated in ISGWM. The actual pumping rates are estimated based on the crop irrigation requirement and well properties, such as the efficiency, pumping capacity and irrigated area. The mean monthly actual irrigation pumping volumes by month are plotted in Figure 4.13. Table 4.4 shows the parameter values for the first 10 wells of the largest pumping rates. In this simulation, the wells were constructed at different times but they all remain active until the end of the simulation. Through the simulation results of ISGWM, it is
135
Annual water depth, mm/a
1000
750
Type 500
a
Precip.
a
Recharge
23.36 % 21.1 %
250
18.37 %
10.24 % 4.62 %
14.15 %
6.24 % 2.42 %
10.12 %7.67 % 6.42 %6.76 %6.39 %
9.15 %
0 2000
2005
Year
Figure 4.11: Annual net recharge versus precipitation. The ratio of the annual recharge and precipitation is labeled along the annual recharge line. found that the crop irrigation requirement estimated according to soil water deficit is significantly larger than those by crop water stress, because the majority of soil water consumption is for soil evaporation. Well 28, although having the largest pumpage, irrigates an area of 0.21 km2 with a pumping capacity of 3815.7 m3 /d which is smaller than other wells. It reveals that the HRU to which Well 28 is assigned undergoes more severe water stress than other areas and it requires more water management practice. Figure 4.13 shows the seasonal changes of irrigation needs. The pattern is controlled by crop growth and water availability. In SWAT, there is a crop growth model to simulate the effects of water availability on crops. Consequently, the crop-water interaction is also simulated in ISGWM. Because the corn is the major crop type in this area, the estimated irrigation water use corresponds well to its
136
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
4572000
4568000
4564000
4560000
4572000
Latitude, m
Recharge mm/day 4568000 1 4564000 0 4560000
−1 Sep
Oct
Nov
Dec
4572000
4568000
4564000
4560000
710000 714000 718000 710000 714000 718000 710000 714000 718000 710000 714000 718000
Longitude, m
Figure 4.12: Mean monthly net recharge distribution. growth cycle in this area as shown in Kranz et al. (2008).
137 Table 4.4: Properties of the first 10 wells with largest pumping rates Starting day
Ending day
Deficit method
Cap.
Eff.
m3 /d
Mean monthly pumping volume, m3
-1843 591 1110 -6440 -4644 -2665 -6348 -13469 -6014 2174
4383 4383 4383 4383 4383 4383 4383 4383 4383 4383
Soil Soil Soil Soil Soil Soil Soil Soil Soil Soil
3815.7 6541.2 4633.3 4088.2 4905.9 4905.9 4633.3 6541.2 5178.4 4088.2
0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
Irrigated area km2
Relative AWC
Irrigated time day
Actual pumping m3 /d
0.21 0.63 0.35 0.69 0.57 0.57 0.38 0.53 0.32 0.38
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
679.4 593.7 471.7 391.1 328.3 328.3 269.7 259.1 229.6 209.0
80000
60000
40000
20000
0 1
2
3
4
5
6
7
8
9
10
11
12
Month
Figure 4.13: Mean monthly irrigation pumping volume for different months.
4.5 Discussion and Conclusion 4.5.1
Innovations
Traditional land surface and hydrologic models oversimplify subsurface water processes and ignore the lateral groundwater flow. Groundwater flow models, on the other hand, cannot depict the physical mechanisms of land surface processes which
138 are often simplified as a constant parameter and need to be externally specified. In this study, an integrated surface-ground water model is developed by coupling SWAT and MODFLOW using a soil water module (SWM). SWM is a non-iterative, mass-conservative and auto-stepping solution to the one dimensional Richards equation. SWM explicitly represents infiltration, soil evaporation, unsaturated water flow, root water update, and lateral drainage. Water-exchanges across the surface-subsurface and unsaturated-saturated zone interfaces are simulated as the system’s dependent top and bottom boundaries of the soil profile, respectively. Stream-aquifer interaction is simulated by the SFR package in MODFLOW. Taking advantage of the simulation capacities of these models, ISGWM can simulate the physical hydrologic processes in three domains and their interactions. The model can also quantify different components in the hydrologic cycle, such as evaporation, transpiration, and groundwater recharge. Irrigation requirement is also simulated using the soil water deficit or crop water stress. By coupling SWAT and MODFLOW, ISGWM inherits the merits of SWAT and MODFLOW. On the other hand, it overcomes some of limitations in these two models. The advantages of ISGWM include: 1. ISGWM estimates the groundwater net recharge considering the effects of land surface processes, soil water movement and the fluctuation of groundwater table. It gives a comprehensive physical description of the water movement during the infiltration and percolation processes. 2. Compared with other integrated three dimensional model, ISGWM is not only computationally efficient, but also is supported with straightforward data pre- and post- processing tools. 3. Actual irrigation pumping estimation accounts for the irrigation deficit and
139 well properties such as well status, pumping capacity and irrigated area. ISGWM is a comprehensive and efficient modeling tool for integrated water resources management. It greatly improves the representation of the unsaturated zone which had been long overlooked or oversimplified.
4.5.2
Future developments
ISGWM has been successfully applied to the Johnson Creek watershed, which is located within the Lower Platte Basin. The recharge process and irrigation activity are simulated both spatially and temporally. It produces reasonable streamflow, recharge and irrigation results. The model can also be applied to evaluation of human activity, land use and climate change effects on not only the surface water but also the groundwater resources. It can be used as a modeling tool for integrated water resources management. However, there are some possible improvements that can be implemented in the future: (1) account for the surface water irrigation and the priorities of the water uses; (2) implement remotely sensed data, such as LAI or ET, in the model
140
Appendix A STREAM DEPLETION RESULTS OF THE LPR MODEL Plate River Depletion with Pumping at ROW 18 COL 95
Plate River Depletion with Pumping at ROW 19 COL 95
Method
Method
SDA TRM
TRM
60
SDF, %
SDF, %
60
SDA
40
20
40
20
Mean of last five year SDA = 23.25% Mean diff. = 4.37%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 26.07% Mean diff. = 3.98%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 20 COL 105
40
TRM
TRM
15
SDF, %
SDF, %
SDA
20
10
10
5
Mean of last five year SDA = 18.93% Mean diff. = 0.55%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 12.29% Mean diff. = −0.51%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 20 COL 75
Plate River Depletion with Pumping at ROW 20 COL 85
120
Method
Method
SDA
SDA
TRM
90
SDF, %
SDF, %
2000
Method
SDA
90
1990
Plate River Depletion with Pumping at ROW 20 COL 110
20
Method
30
1980
Year
60
30
TRM
60
30
Mean of last five year SDA = 27.20% Mean diff. = 6.90%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 27.77% Mean diff. = 7.68%
0 1950
1960
1970
1980
Year
1990
2000
141
Plate River Depletion with Pumping at ROW 20 COL 95 80
Plate River Depletion with Pumping at ROW 22 COL 110
Method
Method
SDA
20
TRM
SDA TRM
60
SDF, %
SDF, %
15 40
20
10
5
Mean of last five year SDA = 28.52% Mean diff. = 3.64%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 13.51% Mean diff. = 0.18%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 22 COL 85
Plate River Depletion with Pumping at ROW 23 COL 110
Method
Method
SDA
20
TRM
SDA TRM
100
SDF, %
SDF, %
15
10
50 5
Mean of last five year SDA = 34.41% Mean diff. = 7.86%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 13.68% Mean diff. = 0.47%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 23 COL 60
Plate River Depletion with Pumping at ROW 25 COL 105
Method
Method 60
SDA TRM
SDA TRM
SDF, %
SDF, %
40 40
20 20
Mean of last five year SDA = 19.93% Mean diff. = 5.32%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 34.21% Mean diff. = 0.99%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 25 COL 120 0.3
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 25 COL 130
Method
Method 0.4
SDA TRM
SDA TRM
SDF, %
SDF, %
0.3 0.2
0.2
0.1 0.1
Mean of last five year SDA = 0.20% Mean diff. = −0.01%
0.0 1950
1960
1970
1980
1990
Mean of last five year SDA = 0.28% Mean diff. = −0.02%
0.0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 25 COL 140 SDA
2000
Method SDA
40
TRM
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 25 COL 60 50
Method 0.3
1980
Year
0.2
TRM
30
20
0.1 10
Mean of last five year SDA = 0.24% Mean diff. = −0.02%
0.0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 20.51% Mean diff. = 4.48%
0 1950
1960
1970
1980
Year
1990
2000
142
Plate River Depletion with Pumping at ROW 25 COL 75
Plate River Depletion with Pumping at ROW 25 COL 85
Method
Method
SDA
SDA 150
TRM
TRM
SDF, %
SDF, %
100
100
50 50
Mean of last five year SDA = 37.51% Mean diff. = 5.47%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 59.33% Mean diff. = 5.31%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 25 COL 90
1990
2000
Plate River Depletion with Pumping at ROW 25 COL 95
Method
Method
SDA 150
1980
Year
SDA
TRM
TRM
SDF, %
SDF, %
100
100
50 50
Mean of last five year SDA = 63.49% Mean diff. = 3.46%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 52.24% Mean diff. = 2.87%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 26 COL 120
Plate River Depletion with Pumping at ROW 26 COL 130 10.0
Method
Method
SDA
SDA
TRM
TRM
7.5
SDF, %
SDF, %
9
6
3
5.0
2.5
Mean of last five year SDA = 7.00% Mean diff. = −0.58%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 5.37% Mean diff. = −0.47%
0.0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 27 COL 120
Plate River Depletion with Pumping at ROW 27 COL 130
Method
Method
SDA
SDA
TRM
TRM
20
SDF, %
SDF, %
10
10
5
Mean of last five year SDA = 22.38% Mean diff. = −1.29%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 10.44% Mean diff. = −0.98%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 27 COL 140 6
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 28 COL 120
Method
Method 40
SDA TRM
SDA TRM
30
SDF, %
SDF, %
4
20
2 10
Mean of last five year SDA = 3.58% Mean diff. = −0.17%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 35.91% Mean diff. = −1.62%
0 1950
1960
1970
1980
Year
1990
2000
143
Plate River Depletion with Pumping at ROW 28 COL 130
Plate River Depletion with Pumping at ROW 28 COL 140 10.0
Method
Method
SDA 15
SDA
TRM
TRM
SDF, %
SDF, %
7.5
10
5
5.0
2.5
Mean of last five year SDA = 15.90% Mean diff. = −1.59%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 7.33% Mean diff. = −0.36%
0.0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 30 COL 105 80
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 30 COL 110
Method
Method
SDA
SDA
60
TRM
TRM
SDF, %
SDF, %
60
40
40
20 20
Mean of last five year SDA = 60.85% Mean diff. = −0.32%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 58.59% Mean diff. = −0.42%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 30 COL 120
1990
2000
Plate River Depletion with Pumping at ROW 30 COL 130
Method 60
1980
Year
Method
SDA
SDA
TRM
TRM
40
SDF, %
SDF, %
20
10 20
Mean of last five year SDA = 58.14% Mean diff. = −1.89%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 25.59% Mean diff. = −2.62%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 30 COL 140
Plate River Depletion with Pumping at ROW 30 COL 50
20 Method
Method
15
SDA 15
SDA
TRM
TRM
SDF, %
SDF, %
10 10
5 5
Mean of last five year SDA = 15.60% Mean diff. = −0.89%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 7.26% Mean diff. = 2.30%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 30 COL 60
Plate River Depletion with Pumping at ROW 30 COL 90 150
Method 40
Method
SDA
SDA
TRM
TRM 100
SDF, %
SDF, %
30
20
50 10
Mean of last five year SDA = 25.13% Mean diff. = 2.29%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 88.23% Mean diff. = 1.26%
0 1950
1960
1970
1980
Year
1990
2000
144
Plate River Depletion with Pumping at ROW 35 COL 50
Plate River Depletion with Pumping at ROW 45 COL 10
Method
Method
SDA TRM
TRM
0.9
SDF, %
SDF, %
15
SDA
10
5
0.6
0.3
Mean of last five year SDA = 7.94% Mean diff. = 1.11%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 0.48% Mean diff. = 0.00%
0.0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 45 COL 75 30
SDA
SDF, %
SDF, %
TRM 20
40
10
20
Mean of last five year SDA = 42.07% Mean diff. = −3.27%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 19.54% Mean diff. = −2.83%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 49 COL 10
1990
2000
Plate River Depletion with Pumping at ROW 49 COL 75 Method
SDA
SDA
TRM
TRM
15
SDF, %
SDF, %
1980
Year
Method
5.0
10
2.5
5
Mean of last five year SDA = 4.62% Mean diff. = −0.08%
0.0 1950
1960
1970
1980
1990
Mean of last five year SDA = 11.78% Mean diff. = −1.95%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 50 COL 20 20
2000
Method
SDA TRM
7.5
1990
Plate River Depletion with Pumping at ROW 48 COL 75
Method 60
1980
Year
Plate River Depletion with Pumping at ROW 50 COL 30
Method
15
SDA
Method SDA
TRM
TRM
15
SDF, %
SDF, %
10 10
5 5
Mean of last five year SDA = 7.91% Mean diff. = 1.16%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 6.27% Mean diff. = −0.30%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 50 COL 40
15
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 50 COL 75
Method
Method
SDA
SDA
TRM
TRM 4
SDF, %
SDF, %
10
2
5
Mean of last five year SDA = 5.97% Mean diff. = 0.06%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 2.30% Mean diff. = −0.40%
0 1950
1960
1970
1980
Year
1990
2000
145
Plate River Depletion with Pumping at ROW 50 COL 85
Plate River Depletion with Pumping at ROW 52 COL 172
100 Method
Method
SDA TRM
TRM
9
SDF, %
75
SDF, %
SDA
50
6
3
25
Mean of last five year SDA = 82.08% Mean diff. = 2.14%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 9.47% Mean diff. = −0.80%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 53 COL 40
1990
2000
Plate River Depletion with Pumping at ROW 55 COL 40
Method
Method 20
SDA 15
1980
Year
SDA
TRM
TRM
SDF, %
SDF, %
15 10
10
5
5
Mean of last five year SDA = 7.80% Mean diff. = 0.03%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 9.29% Mean diff. = 0.00%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 56 COL 170
TRM
75
SDF, %
SDF, %
SDA
TRM
10
50
5
25
Mean of last five year SDA = 15.72% Mean diff. = −1.39%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 84.58% Mean diff. = 3.82%
0
2000
1950
1960
1970
Year
1990
2000
Plate River Depletion with Pumping at ROW 98 COL 105
Method
Method 15
SDA
SDA
TRM
TRM
SDF, %
SDF, %
1980
Year
Plate River Depletion with Pumping at ROW 70 COL 85
40
20
10
5
10
Mean of last five year SDA = 31.00% Mean diff. = 1.49%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 12.52% Mean diff. = 1.09%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 98 COL 140
Plate River Depletion with Pumping at ROW 99 COL 105
Method
8
SDA 7.5
2000
Method
SDA
30
1990
Plate River Depletion with Pumping at ROW 70 COL 130
100
Method 15
1980
Year
TRM
Method SDA TRM
SDF, %
SDF, %
6
5.0
2.5
4
2
Mean of last five year SDA = 5.46% Mean diff. = 0.93%
0.0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 5.41% Mean diff. = 0.44%
0 1950
1960
1970
1980
Year
1990
2000
146
Plate River Depletion with Pumping at ROW 100 COL 10
Plate River Depletion with Pumping at ROW 100 COL 105
Method
Method
SDA 15
SDA
4
TRM
TRM
SDF, %
SDF, %
3 10
2
5 1
Mean of last five year SDA = 10.73% Mean diff. = 0.36%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 1.84% Mean diff. = 0.14%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 100 COL 130
Plate River Depletion with Pumping at ROW 100 COL 140 6
Method
5
Method
SDA
SDA
TRM
4
TRM
SDF, %
SDF, %
4 3
2
2
1
Mean of last five year SDA = 2.20% Mean diff. = 0.29%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 2.69% Mean diff. = 0.58%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 100 COL 3
Plate River Depletion with Pumping at ROW 100 COL 30 40
Method SDA
Method SDA
TRM
TRM 30
SDF, %
SDF, %
40
20
20 10
Mean of last five year SDA = 33.24% Mean diff. = 1.55%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 14.19% Mean diff. = 3.50%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 100 COL 75 SDA
SDA
TRM
SDF, %
SDF, %
TRM
7.5
2
Mean of last five year SDA = 3.04% Mean diff. = 0.34%
0
5.0
2.5
Mean of last five year SDA = 6.28% Mean diff. = 0.52%
0.0 1960
1970
1980
1990
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 100 COL 95
Plate River Depletion with Pumping at ROW 101 COL 10
Method
Method 15
SDA
SDA TRM
SDF, %
TRM
SDF, %
2000
Method
4
15
1990
Plate River Depletion with Pumping at ROW 100 COL 90
10.0
Method
1950
1980
Year
10
10
5
5
Mean of last five year SDA = 11.53% Mean diff. = 2.41%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 9.68% Mean diff. = 0.42%
0 1950
1960
1970
1980
Year
1990
2000
147
Plate River Depletion with Pumping at ROW 101 COL 144 15
Plate River Depletion with Pumping at ROW 103 COL 10
Method
12
SDA
Method SDA
TRM
TRM
SDF, %
8
SDF, %
10
5
4
Mean of last five year SDA = 11.43% Mean diff. = 1.65%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 7.56% Mean diff. = 0.49%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 103 COL 3 30
TRM
SDA TRM
SDF, %
SDF, %
2000
Method
SDA
10
20
10
5
Mean of last five year SDA = 15.56% Mean diff. = 1.10%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 11.80% Mean diff. = 3.21%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 104 COL 142
Plate River Depletion with Pumping at ROW 104 COL 3
Method
Method
SDA 15
1990
Plate River Depletion with Pumping at ROW 103 COL 30
Method 15
1980
Year
SDA
TRM
TRM
SDF, %
SDF, %
10
10
5 5
Mean of last five year SDA = 13.16% Mean diff. = 2.23%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 11.57% Mean diff. = 0.94%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 104 COL 30 30
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 105 COL 10
Method
10.0
Method
SDA
SDA
TRM
TRM 7.5
SDF, %
SDF, %
20 5.0
10 2.5
Mean of last five year SDA = 10.35% Mean diff. = 2.91%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 6.28% Mean diff. = 0.52%
0.0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 105 COL 141
2000
Method
SDA
SDA
TRM
7.5
10
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 105 COL 3
Method 15
1980
Year
5
TRM
5.0
2.5
Mean of last five year SDA = 12.22% Mean diff. = 2.01%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 8.00% Mean diff. = 0.75%
0.0 1950
1960
1970
1980
Year
1990
2000
148
Plate River Depletion with Pumping at ROW 105 COL 30
Plate River Depletion with Pumping at ROW 105 COL 50
Method
Method
SDA
SDA
TRM
TRM
10
SDF, %
SDF, %
20
5
10
Mean of last five year SDA = 7.89% Mean diff. = −1.25%
0
Mean of last five year SDA = 9.15% Mean diff. = 2.65%
0 1950
1960
1970
1980
1990
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 105 COL 95
Plate River Depletion with Pumping at ROW 106 COL 95 5
6 Method
Method
SDA
SDA 4
TRM
TRM
SDF, %
SDF, %
4 3
2
2 1
Mean of last five year SDA = 2.77% Mean diff. = 0.43%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 2.05% Mean diff. = 0.27%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 107 COL 3
2000
Method
SDA
SDA
TRM
TRM
10
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 107 COL 50
Method
2
1980
Year
1
5
Mean of last five year SDA = 7.47% Mean diff. = −1.27%
0
Mean of last five year SDA = 2.04% Mean diff. = 0.27%
0
−5 1950
1960
1970
1980
1990
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 107 COL 95
Plate River Depletion with Pumping at ROW 110 COL 105 3
4
Method
Method
SDA
SDA
TRM
TRM 2
SDF, %
SDF, %
3
2
1 1
Mean of last five year SDA = 1.49% Mean diff. = 0.17%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 0.38% Mean diff. = 0.00%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 110 COL 130
Plate River Depletion with Pumping at ROW 110 COL 140
Method
10.0
SDA
Method SDA
TRM
TRM
4
SDF, %
SDF, %
7.5
5.0
2 2.5
Mean of last five year SDA = 2.33% Mean diff. = 0.38%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 6.72% Mean diff. = 1.08%
0.0 1950
1960
1970
1980
Year
1990
2000
149
Plate River Depletion with Pumping at ROW 110 COL 3
2.0
Plate River Depletion with Pumping at ROW 110 COL 30
Method
Method
SDA TRM
1.5
TRM 10
SDF, %
SDF, %
SDA
1.0
5 0.5
Mean of last five year SDA = 1.27% Mean diff. = 0.25%
0.0 1950
1960
1970
1980
1990
Mean of last five year SDA = 4.09% Mean diff. = 1.38%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 110 COL 95 3
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 111 COL 198 12
Method
Method
SDA
SDA
TRM
TRM
9
SDF, %
SDF, %
2 6
1 3
Mean of last five year SDA = 0.51% Mean diff. = 0.04%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 10.78% Mean diff. = −0.97%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 111 COL 200 5
2000
Method
SDA
SDA
TRM
TRM
6
3
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 112 COL 140
8
Method
4
1980
Year
2
4
2 1
Mean of last five year SDA = 3.65% Mean diff. = −0.32%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 4.29% Mean diff. = 0.77%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 112 COL 196
20
1990
2000
Plate River Depletion with Pumping at ROW 114 COL 143
Method
Method
SDA
SDA 20
TRM
15
1980
Year
TRM
SDF, %
SDF, %
15 10
10
5 5
Mean of last five year SDA = 18.10% Mean diff. = −1.52%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 18.32% Mean diff. = 2.60%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 115 COL 140 8
1990
2000
Plate River Depletion with Pumping at ROW 115 COL 40 40
Method SDA
Method SDA
TRM
TRM
6
30
SDF, %
SDF, %
1980
Year
4
2
20
10
Mean of last five year SDA = 4.48% Mean diff. = 0.81%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 12.72% Mean diff. = 4.71%
0 1950
1960
1970
1980
Year
1990
2000
150
Plate River Depletion with Pumping at ROW 116 COL 139 8
Plate River Depletion with Pumping at ROW 117 COL 138 8
Method
Method
SDA
SDA
TRM
TRM 6
SDF, %
SDF, %
6
4
2
4
2
Mean of last five year SDA = 4.44% Mean diff. = 0.81%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 4.43% Mean diff. = 0.82%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 124 COL 144 20
10.0
2000
Method SDA
TRM
TRM
15
7.5
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 125 COL 140
Method SDA
10
5.0
5
2.5
Mean of last five year SDA = 14.86% Mean diff. = 3.16%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 7.04% Mean diff. = 1.60%
0.0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 126 COL 139
Plate River Depletion with Pumping at ROW 127 COL 138
Method
8
1980
Year
Method
6
SDA
SDA
TRM
TRM
6
SDF, %
SDF, %
4 4
2 2
Mean of last five year SDA = 5.39% Mean diff. = 1.24%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 3.70% Mean diff. = 0.87%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 140 COL 155
1990
2000
Plate River Depletion with Pumping at ROW 140 COL 158
15
8 Method
Method
SDA
SDA
TRM
6
1980
Year
TRM
SDF, %
SDF, %
10
4
5 2
Mean of last five year SDA = 4.88% Mean diff. = 0.11%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 10.41% Mean diff. = 0.08%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 140 COL 160 20
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 150 COL 160
Method
Method 6
SDA TRM
SDA TRM
SDF, %
SDF, %
15
10
4
2 5
Mean of last five year SDA = 15.50% Mean diff. = −0.02%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 3.62% Mean diff. = −0.01%
0 1950
1960
1970
1980
Year
1990
2000
151
Plate River Depletion with Pumping at ROW 150 COL 165
12
Plate River Depletion with Pumping at ROW 150 COL 166
Method
15
Method
SDA
SDA
TRM
9
TRM
SDF, %
SDF, %
10 6
5 3
Mean of last five year SDA = 8.24% Mean diff. = −0.18%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 11.90% Mean diff. = −0.24%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 150 COL 167 20
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 150 COL 170
Method
Method
SDA
SDA 30
TRM
TRM
SDF, %
SDF, %
15
10
20
10
5
Mean of last five year SDA = 15.75% Mean diff. = −0.30%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 27.84% Mean diff. = −0.50%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 160 COL 160 6
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 160 COL 165
Method
Method
6
SDA
SDA
TRM
TRM
4
SDF, %
SDF, %
4
2
2
Mean of last five year SDA = 3.00% Mean diff. = 0.08%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 3.43% Mean diff. = −0.02%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 160 COL 168
Plate River Depletion with Pumping at ROW 160 COL 170
15 Method
Method 20
SDA
SDA
TRM
TRM 15
SDF, %
SDF, %
10
10
5 5
Mean of last five year SDA = 11.46% Mean diff. = −0.20%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 18.61% Mean diff. = −0.39%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 160 COL 180
1990
2000
Plate River Depletion with Pumping at ROW 170 COL 170
Method 60
1980
Year
Method
SDA
SDA
TRM
TRM
40
SDF, %
SDF, %
4
2 20
Mean of last five year SDA = 60.33% Mean diff. = −0.87%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 2.88% Mean diff. = 0.01%
0 1950
1960
1970
1980
Year
1990
2000
152
Plate River Depletion with Pumping at ROW 170 COL 173
Plate River Depletion with Pumping at ROW 170 COL 174
8
12 Method
Method
SDA TRM
9
SDF, %
6
SDF, %
SDA
4
2
TRM
6
3
Mean of last five year SDA = 5.05% Mean diff. = 0.00%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 8.96% Mean diff. = 0.00%
0
2000
1950
1960
1970
Year
Plate River Depletion with Pumping at ROW 170 COL 175 15
2000
Method
SDA
SDA 30
10
SDF, %
SDF, %
1990
Plate River Depletion with Pumping at ROW 170 COL 180
Method TRM
5
TRM
20
10
Mean of last five year SDA = 12.93% Mean diff. = −0.01%
0 1950
1960
1970
1980
1990
Mean of last five year SDA = 32.76% Mean diff. = −0.25%
0
2000
1950
1960
1970
Year
1980
1990
2000
Year
Plate River Depletion with Pumping at ROW 179 COL 194 15
1980
Year
Plate River Depletion with Pumping at ROW 179 COL 195
Method
Method
SDA
20
TRM
SDA TRM
15
SDF, %
SDF, %
10
10
5 5
Mean of last five year SDA = 12.40% Mean diff. = 0.08%
0 1950
1960
1970
1980
Year
1990
2000
Mean of last five year SDA = 20.42% Mean diff. = −0.15%
0 1950
1960
1970
1980
Year
1990
2000
153
Appendix B DATA INPUT INSTRUCTIONS FOR MODFLOW-SDA Input to the Stream Depletion Analysis (SDA) Package is read from the file that has type “SDA” in the name file. All variables are free format. Dataset 1 NSEN IACTIVE ICR ICC ICV IHCOF ISC ICB NBZONEMAX NSEN An integer value that can be specified to be positive or negative. The absolute value of NSEN is equal to the total number of scenarios. If NSEN < 0, SDA runs the baseline and writes the coefficients to disk. If NSEN > 0, SDA reads the coefficients from disk before starting the scenario analysis. If NSEN = 0, SDA is inactivated. IACTIVE File unit number for writing and reading the IBOUND array. ICR File unit number for writing and reading the CR array. ICC File unit number for writing and reading the CC array. ICV File unit number for writing and reading the CV array. IHCOF File unit number for writing and reading the HCOF array. ISC File unit number for writing and reading the SC array. ICB File unit number for writing and reading the head-dependent boundary coefficients. NBZONEMAX An integer value equal to the maximum zones for budget calculation.
154 FOR EACH SCENARIO Dataset 2a ScenName ScenName A label for the scenario run. Dataset 2b ListFile ListFile File name for print information. Dataset 2c SFtype SFname SFtype Solver type to be used, e.g. “PCG”, “SIP”, “DE4”, “PCGN” or “GMG”. SFname File name for the solver input; the format of the solver file is identical with MODFLOW specifications.. Dataset 2d NB, BType1, Btype2, ... BtypeN NB An integer value equal to the number of boundary types for flow response calculation. If NB < 0, all the boundary types will be included in flow response calculation. BType Boundary types for flow response calculation. There are NB types to be specified. BType is a three-letter string which is identical with the Ftype in the MODFLOW name file, e.g. “WEL”, “GHB”, “SFR”, etc. For the storage term, “STO” can be used. Dataset 2e ZONFile ZONFile Zone file name for flow response calculation. The zone file format follows the Zone Budget specification Dataset 2f QTOR, IHED, IWEL, IRCH QTOR A real value used in flow calculation outputs. If the flow smaller than QTOR, it will not be written to the output file. This can help dramatically reduce the size of the flow file because some flows are very small. IHED An integer value used as a flag for writing head changes. If IHED < 0, head changes will be written to an unformatted file. If IHED > 0, head
155 changes will be written to a formatted file with the format specification “(10(1X1PE13.5))”. If IHED = 0, head changes will not be written to disk. IWEL An integer value used as a flag for well simulation in the scenario run. If IWEL > 0, the well package is activated in the scenario run. A ”WEL” file will be needed. The format of the well file follows MODFLOW specification. Please note that pumping rate used for the well package is the accretion flow compared to the baseline condition. IRCH An integer value used as a flag for recharge simulation in the scenario run. If IRCH > 0, the recharge package is activated in the scenario run. A ”RCH” file will be needed. The format of the recharge file follows MODFLOW specification. Please note that recharge rate used for the recharge package is the accretion flow compared to the baseline condition. Dataset 2g [HedFile] HedFile File name for writing head changes. It only reads when IHED > 0. Dataset 2h [WelFile] WelFile Input file name for the well package. It only reads when IWEL > 0. Dataset 2i [RCHfile] RCHfile Input file name for the recharge package. It only reads when IRCH > 0.
156
Appendix C DATA INPUT INSTRUCTIONS FOR ISGWM Input to the unsaturated zone and watershed model can be activated by including a record in the MODFLOW name file using the file type (Ftype) “USW”. The USW file includes specification of the connection of SWAT HRUs and MODFLOW cells, SWAT subbasins and MODFLOW SFR segments, as well as irrigation properties. The format specification of the USW file is explained as follows: Dataset 1 NFRAC NSEG NWEL NFRAC An integer value equal to the total number of interface fractions. NSEG An integer value equal to the total number of SFR segments receiving overland flow from the watershed model. NSEG can be smaller than the total number of SFR segments. In this case, some of the segments do not receive overland runoff but they will still be included in streamflow routing calculation. NWEL An integer value equal to the total number of irrigation wells. FOR EACH FRACTION Dataset 2 HRU ROW COL AREA
157 HRU An integer value equal to the HRU index for the interface fraction. ROW An integer value equal to the row index for the interface fraction. COL An integer value equal to the column index for the interface fraction. AREA A real value equal to the area of the interface fraction. FOR EACH STREAM SEGMENT Dataset 3 SUB SEG SUB An integer value equal to the subbasin index from which the segment will receive overland runoff. SUB An integer value equal to the SFR segment index. FOR EACH IRRIGATION WELL Dataset 4 IH IL IR IC IS IE CAP EFF IDEFI RAWC AREA IRRDT IH An integer value equal to the HRU index in SWAT. IL An integer value equal to the layer index in MODFLOW. IR An integer value equal to the row index in MODFLOW. IC An integer value equal to the column index in MODFLOW. IS An integer value equal to the starting day of the well. IE An integer value equal to the ending day of the well. CAP An real value equal to pumping capacity [L3 /T]. EFF An real value equal to the irrigation efficiency. IDEFI An integer value indicating the method to estimate irrigation deficit. If IDEFI = 1, irrigation deficit is estimated based on transpiration deficit. If IDEFI = 2, irrigation deficit is estimated based on soil water deficit. If IDEFI = 3, irrigation deficit is estimated based on maximum of transpiration or soil water deficit. RAWC An real value equal to the faction of available water content as a threshold for soil water deficit estimation. AREA An real value equal to irrigated area [L2 ]. IRRDT An real value equal to the irrigation time to fulfill the soil water deficit [L].
158
Appendix D R SCRIPTS TO GENERATE SWM AND SUW INPUTS setwd("D:\\@Workspace\\Johnson_Cr\\john_model\\TxtInOut") elev
