UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 1 DEVELOPMENT AND APPLICATION OF NEW ANALYTICAL EXPRESSIONS FOR DRAINABLE AND FILLABLE POROSITY IN SHALLOW UNCONFINED AQUIFERS By S UBODH A CHARYA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2012 PAGE 2 2 2012 Subodh Acharya PAGE 3 3 To my parents, wife and sisters PAGE 4 4 ACKNOWLEDGMENTS Firstly, I would like to thank my parents who instilled in me the very foundations of a moral human being and a decent student. I am also deeply indebted to my sisters Sudha and Sujata for their love and support. I also want to thank all my family members, especia lly my uncle and aunt who supported m e in various ways throughout my student life. A lot of thanks to my cousins Kishor and Anurag, and my friend Sanjay with whom I grew up and have shared the strong est bond of friendship for the past 20 years I thank my advisor Dr. Rao Mylavarapu for his support and guidance throughout my study. I am deeply grateful to him for believing in what I was doing and creating a truly independent working environment for me. I also want to thank my supervisory committee for their valuable suggestions and help to improve my research from its very primordial stage to the final dissertation. Special thank goes to Dr. James Jawitz who gave his valuable time to improve my research I would also like to thank Kelley Hines and Dakshina Mu rthy for their help in sample collection and lab analyses I want to thank the farm crew at the Hastings Research and Education Center for providing all the help I needed with my field work. S pecial thank s to Bart Herrington for his selfless assistance in getting permission from the growers to us e their fields for this research T hank s to all my Nepalese colleagues, seniors and friends at the University of Florida as well as in Gainesville with whom I have enjoyed the compatriotism and the great Nepalese culture. The past six years in Gainesville ha ve been simply amazing. Finally, countless thanks to my beautiful wife Shweta ; without her love and support I would be nowhere near where I am now Her endless love ha s helped me sail through the roughest tides of time and be the best person I can be Thank you my darling. PAGE 5 5 TABLE OF CONTENTS Page ACKNOWLEDGMENTS ................................ ................................ ................................ ............... 4 LIST OF TABLES ................................ ................................ ................................ ........................... 7 LIST OF FIGURE S ................................ ................................ ................................ ......................... 8 ABSTRACT ................................ ................................ ................................ ................................ ... 11 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .................. 13 Unconfined Aquifer Storage Parameters: Drainable and Fillable Porosity ............................ 14 Applications of d and f ................................ ................................ ................................ ......... 18 Modeling Shallow Water Table Dynamics ................................ ................................ ..... 18 Estimation of Groundwater Evapotranspiration ( ETg ) ................................ .................... 20 Hypotheses and Objectives ................................ ................................ ................................ ..... 23 Outline of Dissertation ................................ ................................ ................................ ............ 24 2 AN ALYTICAL EXPRESSIONS FOR DRAINABLE AND FILLABLE POROSITY OF PHREATIC AQUIFERS UNDER VERTICAL FLUXES DUE TO EVAPOTRANSPIRATION AND RECHARGE ................................ ................................ ... 26 Background ................................ ................................ ................................ ............................. 26 Theory ................................ ................................ ................................ ................................ ..... 29 Saturated and Unsaturated Soil Water Storage ................................ ............................... 29 Total Storage Deficit ................................ ................................ ................................ ....... 31 General Expressions for Drainable and Fillable Porosity ................................ ............... 31 Steady Flow from/to the Water Table (Evaporation, Root Water Uptake and Infiltration) ................................ ................................ ................................ ................... 32 Evaporation and root water uptake ................................ ................................ ........... 32 Evapotranspiration approximation based on water table elevation .......................... 33 Water table recharge during rainfall ................................ ................................ ......... 35 d and f during ET or Recharge from the Water Table ................................ .................. 36 Implementation of d and f in a 1D Groundwater Model ................................ ............... 37 Methods ................................ ................................ ................................ ................................ .. 39 Study Site ................................ ................................ ................................ ......................... 39 Water Table and Weather Data ................................ ................................ ....................... 4 0 Parameter Estimation ................................ ................................ ................................ ....... 40 Boundary Conditions and Numerical Simulation ................................ ............................ 41 Results ................................ ................................ ................................ ................................ ..... 42 Drainable Porosity and Fillable Porosity under ET and Recharge ................................ .. 42 Water Table Dynamics ................................ ................................ ................................ .... 43 Discussions ................................ ................................ ................................ ............................. 45 PAGE 6 6 Conclusions ................................ ................................ ................................ ............................. 46 3 ESTIMATION OF EVAPOTRANSPIRATION FROM DIURNAL WATER TABLE FLUCTUATIONS USING DRAINABLE AND FILLABLE POROSITY ........................... 61 Background ................................ ................................ ................................ ............................. 61 Theory ................................ ................................ ................................ ................................ ..... 67 Estimating ETg from DWTF using d and f ................................ ................................ ... 67 Estimation of d and f ................................ ................................ ................................ ..... 68 Hourly Water Table Slopes ................................ ................................ ..................... 69 ETg Calculations ................................ ................................ ................................ ............. 72 Application of the Modified Method ................................ ................................ ...................... 73 Site Description and Water Table Data ................................ ................................ ........... 73 Estimation of Recharge ................................ ................................ ................................ ... 75 Estimation of Subsuface Flux ( q ) ................................ ................................ .................... 76 Penman Monteith ET ................................ ................................ ................................ ...... 77 Results an d Discussions ................................ ................................ ................................ .......... 78 Hourly ETg Estimates ................................ ................................ ................................ ...... 78 Daily ETg Estimates ................................ ................................ ................................ ........ 80 Conclusi ons ................................ ................................ ................................ ............................. 83 4 MODELING SHALLOW WATER TABLE DYNAMICS UNDER A SUBSURFACE IRRIGATION AND DRAINAGE SYSTEM IN FLORIDA ................................ ................. 92 Background ................................ ................................ ................................ ............................. 92 Materials and Methods ................................ ................................ ................................ ........... 95 Study Site ................................ ................................ ................................ ......................... 95 Model Development ................................ ................................ ................................ ........ 96 Conceptual models: irrigation and drainage ................................ ............................. 96 Governing equations ................................ ................................ ................................ 97 Initial and boundary conditions ................................ ................................ ................ 99 Numerical solution ................................ ................................ ................................ 100 Water Ta ble Data and Soil Parameters ................................ ................................ .......... 101 Results and Discussions ................................ ................................ ................................ ........ 102 Model Application ................................ ................................ ................................ ......... 102 Model Eval uation and Sensitivity Analysis ................................ ................................ .. 103 Conclusions ................................ ................................ ................................ ........................... 106 5 SUMMARY AND CONCLUSIONS ................................ ................................ ................... 115 APPENDIX: THE WATER TABLE DYNAMICS SIMULATION MODEL: R CODE .......... 118 LIST OF REFERENCES ................................ ................................ ................................ ............. 127 BIOGRAPHICAL SKETCH ................................ ................................ ................................ ....... 137 PAGE 7 7 LIST OF TABLES Table P age 2 1 Original and Modified van Genuchten (VG) soil moisture retention model parameters and alpha ( a G ) for Clay, Loam and Sand, and Ellzey fine sand series.. ................................ ................................ ................................ ................................ 49 2 2 Nash Sutcliffe Coefficient of Efficiency (C eff ) and Root Mean Square Error ( RMSE ) of the model implemented with hydrostatic (common d and f ) and dynamic (different d and f ) expressions ................................ ................................ ......................... 49 3 1 Root Mean Square Error (RMSE) of hourly and daily ETg estimations using single, hydrostatic d and distict, dynamic d and f ................................ ................................ 85 4 1 Comparison measures: Root Mean Square Error (RMSE) and Nash Sutcliffe Coefficient of Efficiency of water table dynamics simulation for the three wells in site 1 in 2010 and 2011. ................................ ................................ ................................ ... 108 PAGE 8 8 LIST OF FIGURES Figure P age 2 1 Conceptual representation of the water table profile and water balance components in a typical unconfined aquifer with shallow phreatic surface, and bounded by a ditch or stream. ................................ ................................ ................................ ................................ 50 2 2 Conceptual representations of change storage deficit ( D s ) due to a unit rise in WT elevation (left), and changes in unsaturated storage ( U s ) and saturate d storage ( S s ) due to a unit decline in WT(right) in an unconfined aquifer system. ................................ 50 2 3 Relationship between ET rate and WT dep th for clay, loam and sand soils under grass type vegetation, based on equations 2 17 and 2 18.. ................................ ................ 51 2 4 Suction at the soil surface ( as a function of WT depth calculated from HYDRUS 1D and the analytical expressions (equations 2 16 and 2 21) .......................... 51 2 5 Study site in northeast Florida showing agricultural land use and typical water table management system with ditches and irrigation furrows. ................................ ................. 52 2 6 Original and modified van Genuchten moisture retention models fitted to the measured date for Ellzey fine sand series. ................................ ...................... 53 2 7 Drainable (left) and fillable (right) porosity as the function of WT depth for clay, sand and loam soils under grass type vegetation during hydrostatic a nd different steady ET rates from WT table based on equations 2 23 and 2 24.. ................................ 54 2 8 Drainable (left) and fillable (right) poros ity of the Ellzey fine sand series soil under potato at different ET demands. ................................ ................................ ......................... 55 2 9 Drainable (left) and fillable (right) porosity of the Ellzey fine sand series soil under potato for different rainfall (P) rates. ................................ ................................ ................. 56 2 10 Simulated and observed hourly water table dynamics in the Well 1 located at the center of the study field in 2010 (top) and 2011 (bottom) spring potato seasons. ............. 57 2 11 Plot showing the relationship between recharge rate and fillable porosity for Ellzey fine sand. The relationship is presented in terms log( f ) against to enhance the visual pattern in the curves. ................................ ................................ ........... 58 2 12 Observed water table rise at the study site during four rainfall events (3/28, 4/25, 5/04, and 5/17) in 2010, compared to predictions from equations 2 30 and 2 31, neglecting groundwater outflow to the ditches. ................................ ................................ 59 PAGE 9 9 2 13 Water table rise predictions during five hypothetical rainfall events in Ellzey soil with initial WT = 120 cm. Lines are from equati ons 2 30 and 2 31 using recharge dependent f ................................ ................................ ................................ ...................... 60 3 1 Diurnal fluctuations in the solar radiation (red) and observed water table depth (blue) at a northeast Florida site. ................................ ................................ ................................ .. 86 3 2 Drainable and fillable porosity of Ellzey sand series under potato in northeast Florida under different potential ET and precipitation rates. ................................ ......................... 86 3 3 Illustration of contribution of recharge induced sudden WT rise during a rainfall event and extraction of the exact change from the observed by subtracting the contribution of recharge ................................ ................................ ....................... 87 3 4 Calculated hourly actual change in wate r table elevation from the observed data after removing the contribution of recharge during four rainfall events in spring 2010. ................................ ................................ ................................ ........................ 87 3 5 Land use of the three counties in northeast Florida and the two field sites with well positions, chosen for the study. ................................ ................................ .......................... 88 3 6 Conceptual representation of the watertable profiles and the flow towards the ditches in the field. ................................ ................................ ................................ ......................... 88 3 7 Time series of the hourly rainfall, Penman Monteith and groundwater ET predicted by using a single drainable porosity (middle plots) and distinct drainable and fillable porosity (bottom plots) from the two field site s ................................ ................................ 89 3 8 Hourly ETg from the HSM method (black dashed line), DM method (red solid line), and PM ET (blue dots) during clear periods (top) and two rainfall events in 2010. ......... 90 3 9 Same as Figure 8, during 2011 in field site 2. ................................ ................................ ... 90 3 10 Scatter plots of daily Penman Monteith ET and predicted ETg from field site 1 in 2010 (left) and field site 2 in 2011 (right) during spring potato seasons. .......................... 91 3 11 Plots showing the variations in hydrostatic d (open circles) and dynamic d (filled diamonds) during WT drawdown after rainfall (top) and dirunal fluctuations in the WT elevation in a field site in 2010. ................................ ................................ .................. 91 4 1 Land use in the St. Johns, Putnam and Flagler counties in northeast Florida and a field site, planted to potato, chosen for the study. ................................ ........................... 109 PAGE 10 10 4 2 Conceptual representation of water table profiles during irrigation phase. Note that only 2 water furrows are depicted here; the actual number of furrows in the fields may be10 30. ................................ ................................ ................................ .................... 110 4 3 Conceptual representation of water table profile during the early drainage phase (top) and late d rainage phase (bottom). The dashed line represents the initial water table elevation after a precipitation event. ................................ ................................ ................ 110 4 4 Observe d (filled circles) and simulated (solid lines) water table rise due to rainfall (vertical bars) and subsequent drawdown (top), and water table rise after onset of irrigation (bottom) in Well 1. ................................ ................................ .......................... 111 4 5 Scatterplots of the observed and simulated hourly water table depths in the three wells in 2010. The red line indicated 1:1 relationship. ................................ .................... 112 4 6 Scatterplots of observed and simulated hourly water table depths in field site 2 in 2011 potato season. The red line indicate s 1:1 relationship. ................................ ........... 112 4 7 Effect of K s on the time required to lower the water table below 50 cm depth from an initial WT depth of 10 cm (left); and time required to raise the WT to 50 cm from an initial depth of 80 cm (right), in a field of 15 ha area (400m x 400 m). .......................... 113 4 8 Effect of furrow depth on the water table drawdown rate during subsurface drainage. .. 113 4 9 Hourly water table depths during irrigation events of an alternative, 50% water saving irrigation regime (12 our irrigation/drainage cycle). ................................ ............ 114 PAGE 11 11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DEVELOPMENT AND APPLICATION OF NEW ANALYTICAL EXPRESSIONS FOR DRAINABLE AND FILLABLE POROSITY IN SHALLOW UNCONFINED AQUIFERS By Subodh A charya August 2012 Chair: Rao Mylavarapu Major: Soil and Water Science In shallow unconfined aquifers, t he response of the phreatic surface to input and output water fluxes is controlled by two distinct storage parameters, drainable ( d ) and fillable porosity ( f ) which are applicable for water table (WT) drawdown and rise respectively. Although it has long been recognized that these parameters are affected by the soil moisture flux in the unsaturated zone, very few studies have attempted incorporating these effects during d and f estimations. Moreover, only the d estimated from hydrostatic soil moisture profile is commonly used in most studies assuming implicitly that the two parameters are equal In this research, separate analytical expressions for d and f we re developed by accounting for the dynamic soil moisture conditions in the unsaturated zone. The new expressions we re then tested by simulating field WT dynamics estimating the groundwater evapotranspiration (ET) from diurnal WT fluctuations and developing a subsurfa ce irrigation and drainage model under a WT management system in Florida It was found that the occurrence of unsaturated zone fluxes due to evapotranspiration ( ET ) and recharge ( R e ) can result in significantly different values of d and f for a given WT elevation even when hysteresis is neglected Implementation of these expressions in the hydraulic groundwater model to simulate the hourly WT dynamics produced significantly improved results as compared to a hydrostatic d The analytical expr essions also PAGE 12 12 enabled the development of a more complete and accurate method to estimate groundwater ET ( ETg ) from diurnal WT fluctuations corroborating the results of simulation The subsurface irrigation and drainage model, implemented with d and f showed that the model was capable of simulating WT drawdown and rise during respective drainage and irrigation events in a crop field Soil moisture in the field is always under a con tinuous redistribution phase, rarely reaching the hydrostatic condition which therefore stresses the critical importance of incorporating the effect s of moisture fluxes in d and f to study the hydrology of shallow phreatic aquifers This research shows that incorporation of soil moisture dynamics in d and f estimations not only improve s our ability to better understand and explain the hydrology but also aid in developing more robust and complete methods to study it. PAGE 13 13 CHAPTER 1 INTRODUCTION Shallow water table (WT) areas commonly exist in many parts of the world where the water table remains at a small depth from the surface during most parts of the year. These areas include wetlands, riparian zones, forests a s wells as intensively cultivated agricultural lands. Regardless of the vegetation, management type, or climate, the hydrology of shallow WT areas is mostly dictated by the dynamics of phreatic surface [ Nachabe, 2002]. The growth and development of plants in these areas are directly influenced by the dynamics of WT [ Bierkens, 1998: Nachabe, 2002] because of its contribution in the soil moisture fluxes due to evaporation and root water uptake [ Nichols, 1993, 1994; Loheide et al., 2005; Naumburg et al., 2005; Shah et al., 2007; Butler et al., 2007; Laio et al, 2009 etc.]. During precipitation the water table can quickly rise to the surface inducing rapid runoff [ Gillham, 1984: Abdul and Gilham, 1989: Troch, 1992; Hilberts et al., 2007 etc.]. In natural environ ments, such rapid changes in WT elevations in response to the different hydrologic fluxes affect the growth and development characteristics of plants thus acting as a force that can control the species composition of the area. In agricultural land settings such shallow WT poses a considerable challenge to crop production practices Since most of the crops cannot sustain high levels of soil moisture, artificial drainage of the land is required for optimum aeration of the root zone ( Ritzema, 1994; Skaggs, 19 91). On the other hand, the nature of the vegetation, characteristics of the aquifers, and human management practices largely affect the behavior of WT itself. It is, therefore, critically important to understand the interactions among the phreatic surface associated hydrologic components, soil, and the plant characteristics that affect WT dynamics in natural as well as managed environments. PAGE 14 14 Unconfined Aquifer Storage Parameters: Drainable and Fillable Porosity Shallow WT fluctuations in unconfined aquife rs are conventionally estimated by means of a storage parameter called the drainable porosity, also termed as specific yield [ Tang and Skaggs, 1975; Nachabe, 2002]. Drainable porosity ( d ) represents the amount of water released by the aquifer when the WT drops by a unit distance [ Bouwer, 1978; Freeze and Cherry, 1979; Neuman 1987; Nachabe, 2002; Hilberts et al., 2005; Laio et al., 2009]. Drainable porosity is the primary storage parameter of hydraulic groundwater models which have been used extensively to study the drainage of agricultural land [e.g. Glover 1964; Moody, 1966; van Schilfgaarde, 1974, 1991; Bouarfa and Zimmer, 2000, etc.] as well as hillslope drainage [ Childs 1971; Brutsaert, 1994; Verhoest et al., 2002; Singh and Jaiswal, 2006; Hilberts e t al., 2005, 2007, etc.]. Analogous to the concept of d another storage parameter of unconfined aquifer can be defined; the fillable porosity ( f ), which is the amount of water that is stored or absorbed by an aquifer per unit rise in water table [ Wils on et al., 1980; Bouwer, 1978; Sophocleous, 1991; Healy and Cook, 2002; Park and Parker, 2008]. Therefore, for every given WT position, both drainable and fillable porosity should be defined because the WT has the potential to move in both directions (upwa rds and downwards ) depending on the conditions at the aquifer boundaries. The se two parameters therefore, respectively dictate recession and rise of the phreatic surface due to evapotranspiration (ET), rainfall and groundwater flow and are integral parame ters of hydraulic groundwater flow models [ e.g., Chapman and Dressler, 1984; Brutsaert, 1994; Hilberts et al., 2005]. It is well accepted in the literature that fillable porosity f is normally smaller than d due to hysteresis in soil moisture retention [ Bouwer and Jackson, 1974; Bouwer, 1978; Sophoclelous, 1991]. However, hysteretic effects have rarely been considered in hydraulic groundwater models treating the two parameters essentially as one [e.g., Troch et al. 2003; Hilberts et al., 2005, Pumo et a l., 2009, etc.]. Stauffer et al. [1992] was the first and perhaps the PAGE 15 15 only study an unconfined groundwater flow model He applied the dependent domain theory of hysteresis [ Mualem, 1984] to a soil mo isture retention curve model and numerically estimate d d and f during model simulation After Stauffer et al [1992], hysteretic storage parameters have seldom been developed or applied in hydraulic groundwater flow models and the d and f have been trea ted as the common, single parameter. It is possibly due to the complex nature of soil moisture hysteresis which is difficult to incorporate in such hydraulic groundwater models. However, under shallow WT, hysteresis effect may not be as significant as comp ared to other factors. It has also been shown by previous studies that hysteretic effects on soil moisture status and flow may not be as significant in such fields with shallow WT and without ponding conditions [ van Dam et al., 1996]. Despite the extensive use, the nature of the application of d in WT dynamics studies, however, has not been consistent Many studies have treated d as a constant parameter [ Glover, 1966; Skaggs, 1980, 1991; Singh and Jaiswal, 2006] although it has long b een recognized that d is a dynamic parameter that depends on the moisture status of the soil profile [ Childs, 1960]. According to Nachabe [2002], d is constant only if the aquifer response is linear, i.e., the volume of drainage changes linearly with WT fluctuation which does not occur in shallow WT environments [ Duke, 1972; Nachabe 2002]. Some studies estimated d by recording the rise in WT during a brief rainfall period and dividing the amount of rainfall by the rise in WT [e.g., Gerla, 1992]. On th e other hand, several studies have also used a variable d estimated using the soil moisture retention curve to study the subsurface lateral flow and WT dynamics [ e.g. Duke, 1972; Bierkens, 1998; Hilberts et al, 2005, 2007; Laio et al., 2009; Pumo et al., 2010, etc.], which is also the most reliable method. PAGE 16 16 It is well known that the WT dynamics in unconfined aquifers is significantly influenced by the capillary fringe [ Cartwright et al., 2006; Hilberts and Troch, 2006; Nielsen and Perrochet, 2000]. Hence the magnitudes of d and f parameters also depend on both the position of the WT and the soil moisture status in the unsaturated zone [ Hillel, 1998]. In an attempt to incorporate these effects, different expressions for d that account for s oil capillary properties have been developed in the past. The earliest of such expressions was given by Duke [1972] who developed an equation for d as a function of WT depth using Brooks Corey (B C) soil water retention model [ Brooks and Corey 1964]. Pan dey et al. [1992] developed an empirical equation to estimate d to apply in subsurface drainage problems. Bierkens [1998] developed another analytical expression for d as a function of WT depth, by integrating the hydrostatic soil moisture profile above the WT. His expression was based on a variant of the van Genuchten [1980] soil water retention curve model. Nachabe [2002] developed an expression for d also using B C retention model. His study suggested that d was significantly affected by soil capill arity and also depended on the duration of drainage. Hilberts et al [2005] developed an expression analogous to equation using the same concepts and applied in a hillslope subsurface storage and WT dynamics model [ Troch et al., 2003; Hilberts et al., 2007]. Despite the difference in their approaches, all of these methods were based on one common assumption: hydrostatic conditions (zero flux) in the unsaturated zone above the water table. This assumption facilitated the analytical integration of the unsaturated soil moisture profile thus the resulting expressions that partially accounted for the effect s of unsaturated zone moisture retention in estimated d values. Although such hydro static moisture profile based expressions have been successfully used in several studies, it is well known that the storage parameters of unconfined aquifers are also significantly affected by moisture flux in the PAGE 17 17 unsaturated zone [ Childs, 1960; Chapman and Dressler, 1984; Chapman, 1995; Tritscher et al., 2000; Brutsae rt, 2005]. Tritscher et al. [2000] used a two dimensional flow model to show that the error associated with estimation of d from static moisture profile may be as high as 35% when the infiltration rates are high. Chapman and Dressler [1984] showed that fillable porosity was directly affected by the vertical recharge rate ( R e ) and approached zero as R e increased. Thus the application of d expressions based on static moisture profile should ideally be limited to hydrostatic conditions. However, this condition rarely holds in most field conditions, unless the soil is very coarse [ Bear and Cheng 2008]. Assumption of hydrostatic moisture profile in shallow WT environments such as wetlands, riparian areas and crop f ields with controlled WT where a significant portion of the ET demand is directly fulfilled by WT [ Loheide et al., 2005; Nachabe, 2005; Butler et al., 2007]. Another important limitation of d expressions based on the hydrostatic assumption is that they re sult in the same value for both d and f for a given WT depth. These parameters have been treated as common in all of the hydraulic groundwater models and their derivatives [e.g., Troch 2003; Hilberts et al., 2005; Pumo et al., 2009, etc.] which may pote ntially lead to erroneous results especially during water table rise scenarios. Recently there have been attempts to incorporate the effect of unsaturated zone moisture fluxes in d of unconfined aquifers. Laio et al [2009] used B C soil moisture retentio n model to develop an approximate analytical expression for d assuming steady steady state soil moisture profile above WT. Their solution was based on steady upward flo w from WT in presence of roots described by an exponential distribution function. The f inal expression was expressed in terms of a critical WT depth term, representing the point of transition from shallow to deep WT. However, this critical WT depth lacked a closed form expression and was approximated by a PAGE 18 18 hypergeometric function. Although th eir method was a step forward towards incorporating the effect of unsaturated zone hydrodynamics in d expressions these authors also failed to investigate the effect of dynamic soil moisture conditions on f Theoretically, f is a parameter related to t he storage deficit of the aquifer, which can potentially behave differently than d in presence of such fluxes because of the inequality in the change in volume of soil moisture during a unit rise and a unit drawdown in WT respectively Drainable porosity describes the volume of water lost when the WT moves down by a unit distance from a given position. On the other hand, fillable porosity is the volume of water that is taken by the aquifer when the WT moves upward from the same given position. Therefore, in presence of unsaturated zone fluxes (ET or R e ) in the soil, these changes will be unequal which can leads to different magnitudes of d and f parameters for the same WT elevation. However, t heoretical investigation of the behavior of aquife r fillable p orosity during such dynamic conditions has never been investigated. Since the shallow WT can fluctuate rapidly even in a short period of time, it is important to consider separate d and f parameters to describe corresponding WT behavior especially in pr esence of unsaturated zone fluxes due to ET and R e Applications of d and f Modeling Shallow Water Table Dynamics The fundamental processes that determine dynamics of WT in aquifers are described by the principles of water movement through the soil. The most rigorous and exact approach to variably saturated porous media [ Pikul et al., 1974; Tang and Skaggs, 1977]. The equation ha s been used extensively to model the saturated unsaturated flow at various spatial scales in natural as well as managed environments [e.g., Freeze, 1971; Feddes et al., 1974; Taylor and Luthin 1969; Khan and Rushton, 1995, 1996; van Dam and Feddes, 2000 etc.]. PAGE 19 19 However, due to high nonlinearity and associated difficulties of solving this equation, simpler approaches are usually sought [ Pikul et al., 1974; Hilberts et al., 2005]. One such assumption most commonly employed in shallow unconfined aquifers is based on the Dupuit Forc heimmer (D F) theory. This theory assumes that in shallow phreatic aquifers, the vertical hydraulic gradient is small enough to be neglected and the flow occurs essentially below the WT in horizontal direction only [ Childs, 1971; van Schilfgaarde et al., 1 974; Skaggs, 1991]. The D F theory holds true when the lateral dimension of the aquifer is very large as compared to its vertical extent which normally results in small inclinations in the water table profile [ van Schilfgaarde, 1974]. When D F assumptions are valid, the lateral soil moisture flux through the soil is governed by the Darcy equation as (1 1) where, q = subsurface lateral flux; h = the height of WT from the reference datum (usually the impermeable barrier at the bottom of the aquifer), K s is the lateral saturated hydraulic conductivity and x is the space coordinate. Conservation of mass requires that Or, (1 2) where, R e ) system. Equation 1 d or f depending on whether the aquifer is in the state of WT drawdown or WT rise respectivel y. During drainage, the direction of flow is towards the ditches, which drives the WT downwards d If the flow is towards the aquifer, = f because the rate of WT rise PAGE 20 20 depends on the storage deficit of the aquifer. Note that the last term in equation 1 2 can either be ET or R e These components have opposite effects on the WT of the aquifer and therefore, they are dictated respectively by d and f However, this particular approach has never been adopted in hydraulic groundwater model ing studies [e.g., Hornberger et al., 1970; Brutsaert, 1994; Skaggs, 1991; Troch et al., 2003 etc.] which may potentially lead to erroneous estimations of WT fluctuations. Estimation of Groundwater Evapotranspiration ( ETg ) In shallow groundwater environme nts such as riparian zone and wetlands, evapotranspiration is a significant component of the water balance [ Dolan et al., 1984; Loheide et al., 2005; Luo and Sophocleous, 2010] where the shallow WT normally supports most of the plant water requirement [ Sh ah et al., 2007; Gribovszki et al., 2007; Lowry and Loheide, 2010]. Consequently a distinct diurnal fluctuation pattern in the water table (WT) hydrograph is usually observed [ Schilling, 2007; Butler et al., 2007; Nachabe et al., 2005; Loheide, 2008]. This diurnal water table fluctuation (DWTF) corresponds to the 24 hour variation cycles in the climatic variables, solar radiation, temperature and humidity [ Gribovszki et al., 2007 ; Gribovszki et al., 2010]. Since there is normally a direct link between the vegetation water uptake and the water table (WT), ET loss from the shallow WT environments may be estimated from the soil water balance. A general expression for ETg from soil water balance can be expressed as (1 3) where, h = water table elevation = the storage parameter of the aquifer q = groundwater flux, PAGE 21 21 R e = rainfall recharge and ETg = groundwater evapotranspiration loss The last term in the right hand side of equation 1 3 gives the total change in soil moisture storage per unit area expressed in term of WT elevation. White [1932] was the first to analyze and report such direct link between ET and WT fluctuations in shallow water table areas. His method, however, was bas ed on different simplifying assumptions. The most critical assumption of original method was that the parameter was assumed constant. He also assumed that the rate of groundwater flow to the area was constant over the 24 hour period and could be estimated by the rate of WT rise between midnight and 4:00 am when ETg is negligible or zero. The original equation for estimating ETg cane be expressed as (1 4) where, d is the drainable porosity (or specific yield) of the aquifer, r is the rate of WT rise between midnight and 4:00 AM when the ET is considered zero, and is the rise or fall of WT during the 24 hour period. Equation 1 4 suggests that, even if the underly ing assumptions hold true, original method is applicable only at daily time scales since it considers constant groundwater flow over the 24 hour period. However, subdaily scale ET estimation may be necessary in many situations which cannot be obtai ned from equation 1 4. Following White many have investigated the DWTF phenomenon to estimate the shallow groundwater consumption by vegetation [e.g., Dolan et al., 1984; Gerla, 1992; Loheide et al., 2005, Nachabe et al., 2005, Lautz, 2008, Butler et al., 2007, Gribovszki et al., 2007; Loheide, 2008; Mould et al., 2010, etc.]. Some of these studies also evaluated the underlying PAGE 22 22 assumptions in equation 1 3 and found that the assumptions of constant groundwater flow as well as constant d are normally not met in field conditions [ Gerla 1992 ; Loheide et al., 2005; Gribovszki et al., 2007]. Gerla [1992] reported that assuming a constant d in equation resulted in excessively large ETg estimation in a wetland environment. According to Loheide et al. [2005], the assumption of constant d creates problems in the White method and showed, with the help of soil moisture flow simulations, that d was critically important in accurately estimating ETg from DWTF. Loheide [2008] proposed a modified White metho d to estimate ETg at subdaily time scales considering variable groundwater flow rates. Gribovszki et al. [2007] presented a similar modification to Loheide [2008] by using successive steady state groundwater flow and a WT dependent d expression developed by Nachabe [2005]. While these studies have shown that implementation of variable groundwater flow and variable d can generally resolve the inherent limitations of the White method and allow for subdaily scale ETg estimations [ Gribovszki et al., 2007; Lo heide, 2008] their study neglected the effect of unsaturated flux in d estimations. As a result, these methods cannot accurately reflect rapid WT fluctuations shown by shallow phreatic aquifers during ET and precipitation. Especially during and after ra infall events, the assumption of hydrostatic moisture profile is severely disrupted because of rapid soil moisture redistribution. Therefore, the actual d value may differ significantly from the estimated hydrostatic d [ Tritscher et al., 2000] consequent ly resulting in potentially erroneous ETg estimations. In order to avoid these complications during rainfall, most of the White based ETg estimation methods presented in the literature typically omitted the days with precipitation in their calculations [e.g. Dolan et al., 1984; Gribovszki et al., 2007; Schilling, 2007]. As a result, these estimations are usually discontinuous at smaller tim e PAGE 23 23 scales (i.e. hourly or daily) although fairly good estimates may be obtained on larger time scales (e.g. monthly average) if the rainfall frequency is low. The potential impact of rainfall on the underlying assumptions of White based ETg estimation meth ods results due to the disturbance in the DWTF signal of the WT. However, the literature provides very little work on the extent to which the ETg estimates are affected when days with rainfall are included in the calculation. In one study, Gerla [1992] inc luded the days with precipitation to estimate ET using the White method. He modified the original White method by accounting for the effect of infiltration on water table rise thus extending the use of the method even during precipitation events. However, h is study not only used a d value estimated from the WT rise during rainfall but also failed to implement separate d and f during ETg calculations. Theoretically, use of these distinct d and f parameters should enable the estimation of ETg from soil moisture balance continuously regardless of the rainfall events because it is possible to estimate the contribution of rainfall to change in WT if reliable estimates of recharge ( R e ) and corresponding f are available. Hypotheses and Objectives The mai n hypotheses of this research are 1. Drainable and fillable porosity of shallow unconfined aquifers can behave in significantly different manner in presence of unsaturated zone vertical fluxes due to recharge and evapotranspiration respectively 2. Dynamics of shallow WT in unconfined aquifers can be understood and modeled significantly better when both drainable and fillable porosity parameters are implemented in the hydraulic groundwater model (Boussinesq equation) 3. Estimates of direct groundwater evapotranspiration based on the White [1932] method can be reliably extended to rainfall periods as well if reliable estimate of fillable porosity is available. The major objectives of the research are PAGE 24 24 1. To develop separate analytical expressions for the drainable and fillable porosity of unconfined aquifer in presence of hydrodynamic conditions and compare with a previous expression estimated under hydrostatic conditions 2. To implement the flux dependent d and f parameters in the hydraulic groundwater model to simulate shallow WT dynamics and compare with the model implemented only with hydrostatic based d parameter 3. To develop a modified White [1932] based method to estimate groundwater ET from diurnal WT fluctuations usin g the flux dependent d and f parameters. 4. To develop and test a shallow WT dynamics model, implementing the flux dependent d and f in order to study the effectiveness of WT control method under a subsurface drainage and irrigation system. Outline of Dis sertation This dissertation presents the development of analytical expressions for the storage parameters d and f and their application in estimation of ETg and development of a WT dyamics model for agricultural fields with WT management system s In C hapter 2, analytical expressions for d and f are developed followin g their p ro per definitions from the literature. Firstly the generalized analytical expressions of the two parameters under hydrodynamic conditions are presented. More specific expressions that account for the steady state unsaturated flow to and from WT are then developed. These expressions are then implemented in the Boussinesq equation to simulate the WT dynamics in a field site in northeast Florida and compared with the field data. Spec ific comparisons are also presented with the Boussinesq equation implemented with the single d parameter estimated from the hydrostatic moisture profile. In C hapter 3, the d and f expressions are applied to estimate ETg using the observed diurnal fl uctuations in WT data. The common White based method of estimating ETg is modified by implementing the two parameters in order to apply the method even during and after rainfall periods which is a significant advantage over the previously available methods PAGE 25 25 Application of the modified method during rainfall periods is enabled by the use of recharge dependent f which allowed for removing the effect of recharge from observed WT fluctuation time series The modified method is then used to estimate ETg from a potato field in northeast Florida during 100 days period over the two year period from 2010 to 2011. The results are then compared with the standard Penman Monteith ET estimates and also the ETg estimates from a previously available single hydrostatic d parameter method. C hapter 4 describes the develop m ent application and evaluation of a WT dynamics model for agricultural fields under a WT management system typically practiced in Florida. Conceptual models of the WT movement during subsurface irrigat ion and drainage via ditches and furrows are presented. The 1D Boussinesq equation that governs the groundwater movement in shallow unconfined aquifer also governs the flow in these ditch drained lands. The governng equation implemented with separate d an d f is solved numerically by implicit finite difference method to simulate WT drawdown during drainage and WT rise during irrigation. The model results are compared to the WT data from two field s during two consecutive potato seasons in 2010 and 2011 Th e sensitivty of the model to some important field and soil parameters is also assessed. Simulated behavior of WT under an alternative potentially more water saving management scenario is also pre sented in this chapter. Finally i n C hapter 5 the summary of the major findings and the main conclusions of the research are presented. PAGE 26 26 CHAPTER 2 ANALYTICAL EXPRESSIONS FOR DRAINABLE AND FILLABLE POROSITY OF PHREATIC AQUIFERS UNDER VERTICAL FLUXES DUE TO EVAPOTRANSPIRATION AND RECHARGE Background In many low lying, humid regions where shallow water tables (WTs) are common, vegetation growth is largely driven by the water table dynamics [B ierkens, 1998; Nachabe 2002]. It is therefore important to understand WT dynamics for efficient management of crops as well as natural vegetation. Shallow WT fluctuations in unconfined aquifers are generally studied by means of a storage parameter called the drainable porosity (or specific yield). Drainable porosi ty ( d ) represents the amount of water released by the aquifer when the water table drops by a unit distance [ Bouwer, 1978; Freeze and Cherry, 1979; Neuman 1987]. It is estimated as the change in total water storage of the aquifer, S per unit drop in WT e levation [ Bear 1972]: (2 1) W here h is height of the phreatic surface above a reference datum. Analogous to the concept of d is the fillable porosity ( f ), which can be defined as the amount of water that is stored or absorbed by an aquifer per unit rise in water table [ Wilson et al., 1980; Bouwer, 1978; Sophocleous, 1991; Park and Parker, 2008]. It can be estimated, in a similar manner to d as the change in total storage deficit of the aquifer, D s per unit rise in WT elevation: (2 2) Drainable and fillable porosity play key roles in WT fluctuations due to evapotranspiration (ET), rainfall, and are integral parameters of Boussinesq type groundwater flow models [ e.g., Chapman and Dressler, 1984; Brutsaert, 1994; Chapman, 2005; Hilberts et al., 2005]. These PAGE 27 27 parameters are also required when ET is estimated from groundwater fluctuations [ Loheide et al., 2005; Gribovszki et al., 2007; Loheide, 2008]. In shallow phreatic aquifers, WT fluctuations are significantly influenced by the soil capillary retention [ Childs, 1960; Gillham, 1984; Nielsen and Perrochet, 2000; Healy and Cook, 2002; Hilberts and Troch, 2006; Cartwright et al., 2006]. Hence the magnitudes of d and f depend on the position of WT and soil moisture status in the unsaturated zone [ Hillel, 1998], rather than being unique for a given soil type [ Hilberts et al., 2005]. While existing analytical expressions for d do account for the effect of unsaturated zone soil moisture [e.g., Duke, 1972; Bierkens, 1998; Nachab e, 2002; Hilberts et al., 2005], these methods assume hydrostatic conditions in the soil profile. Although d expressions derived using static moisture profiles have been successfully used in several studies, it is well known that the storage parameters of unconfined aquifers are also significantly affected by moisture flux in the unsaturated zone [ Childs, 1960; Chapman and Dressler, 1984; Chapman, 1995; Tritscher et al., 2000; Brutsaert, 2005]. Tritscher et al. [2000] used a two dimensional flow model to show that the error associated with estimation of d from static moisture profile may be as high as 35% when t he infiltration rates are high. Chapman and Dressler [1984] showed that fillable porosity was directly affected by the vertical recharge rate ( R e ) and approached zero as the recharge rate increased. Thus the application of d expressions based on static mo isture profile should ideally be limited to hydrostatic conditions. However, this condition rarely holds in most field situations, unless the soil is very coarse [ Bear and Cheng 2008]. The assumption of hydrostatic pressure distribution is especially inap propriate in shallow WT environments such as wetlands, riparian areas and crop fields with controlled WT where a PAGE 28 28 significant portion of the ET demand is directly fulfilled by the water table [ Loheide et al., 2005; Nachabe, 2005; Butler et al., 2007]. Anot her important limitation of d expressions based on the hydrostatic assumption is that they result in the same value for both d and f for a given WT depth. Normally due to hysteresis in moisture retention and release by soils, f is regarded as smaller t han d [ Bouwer and Jackson, 1974; Bouwer, 1978; Sophoclelous, 1991]. As we show below, as the moisture profile deviates from a hydrostatic pressure distribution, different d and f values result for a given WT depth, even when hysteresis is neglected. Yet separate d and f parameters for water table decline and rise have not previously been implemented in unconfined groundwater flow studies. These parameters have also been treated as common in all of the hydraulic groundwater models and their derivatives [e.g., Troch 2003; Hilberts et al., 2005; Pumo et al., 2009, etc.] which may potentially lead to erroneous results especially during water table rise scenarios. Laio et al. [2009] recently developed an approach for incorporating dynamic moisture profile s explicitly into d estimations. These authors used the soil moisture retention model of Brooks and Corey [1964] to develop an approximate analytical expression for d accounting for both recharge and root uptake. Their relation was expressed in terms of a critical WT depth term, representing the point of transition from shallow to deep WT. However, this critical WT depth lacked a closed form expression and was approximated by a hypergeometric function. Also, these authors did not separately consider both d and f The main objective of this study was to develop closed form expressions for d and f as separate aquifer parameters that account for vertical soil moisture flow in the unsaturated zone. Assuming no hysteresis in soil moisture retention, steady state vertical soil moisture fluxes (ET and recharge) were considered in successive times to est imate the moisture profile above the PAGE 29 29 water table. The moisture profile was then integrated to obtain S and D s Here, we use a variant of the van Genuchten [1980] retention model [ Troch 1992; Troch et al. 1993] which enables the analytical solutions for b oth S and D s Differentiating these according to Equations 1 and 2 results in generalized expressions for d and f in terms of the matric suction at the soil surface ( ) and soil hydraulic parameters. This study therefore provides an expression for a s eparate f parameter for a given WT depth, in addition to d which has not been previously investigated. The effects of unsaturated soil moisture flux to and from the water table are accounted for by estimating under successive steady flu x ET or rech arge ( R e ) using the hydraulic conductivity model of Gardner [1958]. The final d and f expressions are flexible enough to incorporate any predefined ET and R e functions. Finally, the analytical expressions are implemented in a one dimensional hydraulic gr oundwater model to simulate shallow WT dynamics in a 15 ha potato field in northeast Florida, managed under a WT control system. Simulated hourly WT dynamics are compared for models implemented with only a single, hydrostatic profile based expression for d and with distinct d and f expressions developed in this study based on dynamic moisture profile. Theory Saturated and Unsaturated Soil Water Storage We consider an unconfined aquifer with a shallow WT and plants growing at the soil surface. The aquifer is assumed to be incompressible with no hysteresis in the soil moisture retention behavior. The water balance components of the cross section of such an aquifer with a shallow WT bounded by a ditch or stream are illustrated in Figure 2 1. A concept ual representation of saturated ( S s ) and unsaturated storage ( U s ), storage deficit ( D s ), and the changes in these quantities due to WT decline or rise is presented in Figure 2 2. The unsaturated soil water storage is normally calculated by integrating the moisture profile from the WT to the soil PAGE 30 30 surface with respect to depth (z). However, analytical integration is possible only when the pressure distribution in the unsaturated zone is hydrostatic so that the matric suction can be expressed as height above t he WT [e.g., Bierkens 1998; Hilberts et al., 2005]. Under dynamic conditions, is a non linear function of z and depends on the magnitude of the flux, precluding analytical integration of the moisture profile with respect to depth. Since soil water conten t in the unsaturated zone is directly related to the soil matric suction ( ) at each respective depth, we integrated with respect to Integrating the moisture profile with respect to from the limit of WT depth ( ) to the soil surface ( ) yields the total unsaturated storage. Therefore, from Figure 2 2 the total soil water storage ( S ) of the profile can be estimated by taking the sum of the area and area under ABDE as (2 3) w here is the matric suction at respective depth z, and are the saturated and residual soil water contents respectively, h is height of the phreatic surface above the reference datum, and corresponds to the matric suction at the soil surf ace. The two terms on the right hand side of equation 2 3 respectively represent the saturated and unsaturated water storage fractions of aquifer vertical cross section. Integration of the second term in equation 2 3 requires a constitutive ( ) relationsh ip. In this study we used a modified van Genuchten [1980] retention model introduced by Troch [1992], which has been used in different studies [e.g., Bierkens, 1998; Hilberts et al., 2005, 2007]: (2 4) where and n are soil specific parameters of the modified model. This modified relationship enables integration of equation 3 which gives PAGE 31 31 (2 5) where and are the parameters of the modified van Genuchten model Total Storage Deficit The total water storage deficit of the aquifer vertical cross section can be estimated by taking the area under BCD in Figure 2 2 as (2 6) Substituting in equation 6 with the modified van Genuchten retention model and integrating yields (2 7) It also follows from equations 2 5 and 2 7 that General Expressions for Drainable and Fillable Porosity From Equations 2 5 and 2 7, d and f can be obtained by taking the derivative of S and D s with respect to h : (2 8) (2 9) It also follows from equations 2 8 and 2 9 that, for a given water table depth, d and f can be related as (2 10) Equations 2 8 and 2 9 are generalized expressions for d and f as functions of the WT elevation. Note that both equations 2 8 and 2 9 can be defined for every WT elevation because at every PAGE 32 32 point below the surface, the phreatic surface has a potential to move either upwards or downwards based on the combined effect of the water balance components (see Figure 2 1). In both equations 2 8 and 2 9, the critical variable is matric suction at the soil surface. Under hydrostatic condition, the magnitude of is simply the height of the soil surface above the WT, i.e., When this is substituted in equations 2 8 and 2 9, we get yielding the drainable porosity expression derived by Bierkens [1998] and Hilberts et al. [2005]: Under dynamic soil moisture conditions, on the other hand, the relationship between and h becomes much more complex. Due to highly nonlinear nature of soil water movement, normally it is not possible to obtain an explicit relationship between h and under transient flow conditions. The simplified approach used here is to assume steady state flux above the water table, which enables direct calculation of Steady Flow from/to the Water Table (Evaporation, Root Water Uptake and Infiltration) In th is section, expressions for are derived considering steady state evaporation or infiltration under shallow water table conditions. Evaporation and root water uptake From Equations 2 5 and 2 7, d and f can be obtained by taking the derivative of S a nd D s with respect to h : (2 8) (2 9) PAGE 33 33 It also follows from equations 2 8 and 2 9 that, for a given water table depth, d and f can be related as (2 10) Equations 2 8 and 2 9 are generalized expressions for d and f as functions of the WT elevation. Note that both equations 2 8 and 2 9 can be defined for every WT elevation because at every point below the surface, the phreatic surface has a potential to move either upwards or downwards based on the combined effect of the water balance components (see Figure 2 1). In both equations 2 8 and 2 9, the critical variable is matric suction at the soil surface. Under hydr ostatic condition, the magnitude of is simply the height of the soil surface above the WT, i.e., When this is substituted in equations 2 8 and 2 9, we get yielding the drainable porosity expression derived by Bierkens [1998] a nd Hilberts et al. [2005]: Under dynamic soil moisture conditions, on the other hand, the relationship between and h becomes much more complex. Due to highly nonlinear nature of soil water movement, normally it is not possible to obtain an explicit relationship between h and under transient flow conditions. The simplified approach used here is to assume steady state flux above the water table, which enables direct calculation of Evapotranspiration approximation based on water table elevation In case of natural shallow WT environments such as wetlands and riparian zones, most of the vegetation is well adapted to perpetual saturated or near saturated conditions, thus showing little or no reduction in root water uptake due to high soil moisture [ Laio et al., 2009; Pumo et PAGE 34 34 al., 2009; Tamea et al., 2009, 2010]. In such environments, it is reasonable to assume that which yields Both linear [ Harbaugh, 2005] and exponential [ Shah et al., 2007] relations have been developed to estimate ET from water table elevation. In this study we used the exponentia l function given by Shah et al. [2007] since it considers both soil and vegetation type: (2 17) (2 18) where and reduction in ET as the WT falls below (Figure 2 1). It was also shown by Shah et al. [2007] that the magnitude of and depend both on the type of soil and vegetation. S maller and larger values correspond to finer soils and deeper rooted vegetation and vice versa. Figure 2 3 illustrates the relationship in 2 17 and 2 18 for clay, loam, and sand soils under typical grass type vegetation based on the parameter values r eported by Shah et al. [2007]. Note that for clay soil the reduction in the rate of ET with decreasing WT is much smaller followed by loam and sand. Estimation of the ET parameters for our field site is discussed in the parameter estimation section. The pa rameter values for the ET vs WT relationship for major soil textural classes under different vegetation types can be found in Shah et al. [2007]. Once the ET is estimated from equations 2 17 and 2 18, can be easily calculated by using equation 2 16. C omparison between the value calculated by equation 2 16 and the values obtained from numerical solution of Richards equation, using HYDRUS 1D model [ Simunek et al., 1998], for estimated ET values for our study area is presented in Figure 2 4. Note that the numerical solution and the PAGE 35 35 steady state analytical solution do not differ significantly, suggesting that our assumption holds nearly true during ET from the water table. Water table recharge during rainfall Water table recharge during rainfall begins once the wetting front reaches the phreatic surface [ Woods et al., 1997] resulting in its rise. Depending on the depth to WT and the soil hydraulic properties, it can be assumed that a fraction of total rainfall during a time step (e.g., daily, hourly) co ntributes to recharge [ Park and Parker, 2008] and subsequent water table rise. Under shallow WT conditions, almost all of the infiltrated water contributes to recharge instantaneously [ Novakowski and Gillham 1988; Hilberts et al., 2007] because the capill ary fringe is usually close to the surface [ Laio et al., 2009]. As the WT depth increases, more infiltrated water is held in the unsaturated zone, thus reducing the fraction of recharge during the same interval of time. A simplified approximation method wo uld be to estimate recharge as a function of precipitation and WT depth without having to explicitly model the unsaturated flow processes. For example, Woods et al. [1997] estimated spatially variable recharge as a function of effective precipitation and WT depth in a hillslope. However, their recharge function was solely dependent on the WT elevation regardless of the soil type. Here, we assume that the recharge is spatially uniform throughout the field site and depends on the soil type, in addition to the position of the phreatic surface. In this study, we propose that the recharge fraction of rainfall at a given time step can be estimated by an exponential function expressed as (2 19) (2 20) where R e is the recharge rate, P is the rainfall amount received during the given time step, and r and d r are the parameters that determine what fraction of rainfall is co nverted to recharge PAGE 36 36 depending on the soil type. Note that the unsaturated flow processes are not explicitly considered when equation 2 19 is used to estimate WT recharge. Therefore, it cannot be effectively applied at very short time steps [Woods, et al., 1997]. However, if the WT is close to the surface, hourly time step may be adequate since instantaneous recharge can be assumed [e.g., Laio et al., 2009] Now, if we treat R e ( ve) as the steady flux in successive times steps, is given as 2 Note that because we make the simplifying assumption of constant flux, R e given by equation 2 21 may not correspond to values measured during rainfall infiltration. is likely to reach zero (or positive) before the water table experiences any recharge. Substituting a positive value of in equations 2 8 and 2 9 would then produce physically impossible negative d or f The values o btained from the HYDRUS 1D model and equations 2 21 are compared in Figure 2 4 for a constant rainfall rate of 0.03 cm/hr at two initial WT depths. The two solutions produced similar values for the shallow water tables when the rainfall rate was small d and f during ET or Recharge from the Water Table Differentiating equation ( 2 16) and ( 2 21) with respect h yields = 2 ET (+ ve) or recharge, R e ( ve). Combining equation 2 22 with equations 2 8 and 2 9 gives the final expressions for d and f under successive steady state fluxes from the water table 2 PAGE 37 37 2 It can be seen from equations 2 23 and 2 24 that d f The disparity between the two parameters is determined by the s nder hydraulic equilibrium where ; and 2 2 2 d and f Implementation of d and f in a 1 D Groundwater Model Movement of water through soils and other porous media are most accurately represented by the Richards equation. However, due to high nonlinearity and associated difficulties of solving this equation, simpler approaches are usually sought [ Pikul et al., 1974]. The overall water balance per unit area of an aquifer can be expressed as (2 25) Note that in equation 2 25 only a single porosity parameter is implemented, which is the conventional approach. In equation 2 25, substitution for the groundwater flux, assuming horizontal flow [ Bear and Cheng, 2008] and neglecting vertical flow, yields the one dimensional Boussinesq eq uation PAGE 38 38 (2 26) where the first term on the right gives lateral the groundwater flux. However, using only a single in equation 2 26 to express both ET and R e in terms of WT elevation is theoretically inconsi stent. Therefore, implementation of separate d and f developed in this study requires reformulation of equation 2 26, depending on the direction of groundwater flow: (2 27) during groundwater inf low towards the aquifer from the lateral boundaries, and (2 28) during groundwater outflow from the aquifer towards the boundaries. Note that only one or both of the storage parameters may appe ar in the governing equation depending on which unsaturated fluxes are active. For example, when groundwater is flowing in to the system and ET = 0, then the governing relation (equation 2 27) will include only fillable porosity. But for ET > 0, both and appear in the governing equation because groundwater inflow raises the WT (controlled by ) while ET tends to lower the WT (controlled by Hence, the response of the phreatic surface will depend on the net difference between the magnitude of groundwater inflow and the ET loss. A common example of this condition is the diurnal fluctuations normally observed in shallow WT environments [e.g., Loheide et al., 2005 ; Gribovszki et al., 2007; Loheide, 2008] where the WT normally drops during the day time (ET > groundwater inflow) and recovers In equation 2 28, note that the recharge component is linked with both and An example scenario governed by equation 28 occurs in artificially drained fields where the water table is usually lowered before anticipated rainfall events in order to avoid the inundation of the PAGE 39 39 field. While the groundwater flow direction at the later al boundaries is outward, the response of the water table (rise or drop) will depend on the magnitude of the recharge at the phreatic surface. When R e > 0, the magnitude of d increases because at each time step there will be some accretion of moisture [ C hapman and Dressler, 1984] which needs to be removed from the aquifer to bring the same decline in the water table as compared to when R e = 0. Therefore, at each time step with R e > 0, the total amount of moisture added by recharge to the system will be pa rtitioned into two fractions: (a) the fraction which contributes to WT rise and is controlled by f and (b) the fraction which is drained from the profile as the WT drops due to groundwater outflow, and is controlled by d Therefore, the net effect of th is R e is determined by the difference between the reciprocals of and Methods Study Site The study field chosen to test the expressions for d and f was located in northeast 5). This area is a low (<10 masl ), flat landscape with a naturally shallow water table (< 2m from surface). A locally important crop is potato grown under a WT management system known as seepage irrigation [ Smajstrla, 2000; Munoz et al., 2008; Acharya and Mylavarapu, 2011] wherein the WT is maintained close to the crop root zone by applying water from shallow furrows spaced uniformly at 18 20m between main lateral ditches (1.5 2.0m deep) at the edges of fields on the order of 3 15 ha. The flow system in these agricultural fields therefore can be appropriately represented by the Boussinesq equation (i.e., equations 2 27 and 2 28). The soil in the field site is highly sandy classified as Ellzey fine sand series (sandy, siliceous, hyperthermic, Arenic Endoaqualf) [ NRCS, 1999] with some clay c ontent at depths 75 120cm [ Acharya and Mylavarapu 2011]. PAGE 40 40 Water Table and Weather Data Three shallow (1.2m deep) wells were installed in a 15 ha field (Figure 2 5) prior to the Spring 2010 potato season (February to June). Hourly WT data were collected from each well using data logging pressure transducers. In the 2011 potato season, hourly WT data were collected from 3 new wells installed in the same field at approximately the same locations. In both years, hourly PET estimates were calculated using the Penman Monteith [ Monteith, 1965] equation based on the weather data collected from a Florida Automated Weather Network [ FAWN, 2011] weather station located adjacent to the field site. Hourly PET estimates were converted to the actual potato ET estimates u sing available potato crop coefficient values [ Allen et al., 1998]. Parameter Estimation Moisture retention parameters of Ellzey fine sand were estimated by fitting soil moisture retention data to the modified van Genuchten function (equation 2 4) which sh owed a very good fit for both the original and modified equations (Figure 2 6). Acharya and Mylavarapu [2011] reported a lab determined average K s = 10.0 cm/hr for Ellzey fine sand, while Rosa [2000] reported K s = 7.0 cm/hr below the WT of the same soil us ing a, slug test method. In this study we used the latter value since in situ observations are more likely to represent field conditions. Table 2 1 shows the values of the hydraulic parameters for Ellzey series. The most rigorous method for estimating parameter would be to fit the model to experimental K data. However, due to the constraints in collecting such experimental data, it is often necessary to revert to other methods. One such method is estimating by establishing correspo ndence with the parameters of the soil moisture retention model [e.g., Zhu et al., 2004; Rucker et al., 2005; Ghezzehei et al., 2007]. Ghezzehei et al. [2007] developed a simple formula to calculate from the original van Genuchten nd n by PAGE 41 41 defining two values at which and match. They reported that, for soils with n > 2, can be estimated reliably by (2 29) Equation 2 29 was therefore used to estimate for Ellzey sand since the value of n for Ellzey series was 2.62 (Table 2 1). Parameters and of the ET WT relationship in 2 17 and 2 18 were adapted from Shah et al. [2007] who determined the parameter values for the major soil textural classes under bare soil, grass and forest type vegetation. The soil in our field site was dominated by fine sand with clay content up to 11% at 70 120cm depth [ Acharya and Mylavarapu, 2011]. Based on these results, the soil was considered close to the loamy sand texture. The parameter values were then determined by taking the average of the values for bare soil and grass since the root system of potato is normally shallow. The estimated values of and for Ellzey sand were 0.08 and 45cm respectively. In order to estimate the parameters of the recharge function given by equation 2 19, several hypothetical simulations of rainfall infiltration with different initial WT depths and rainfal l rates were carried out using HYDRUS 1D. The magnitudes of the soil moisture flux at the saturated unsaturated zone interface (i.e. the phreatic surface) were recorded at each time step which were then then fitted with the recharge function in 2 20, whic h provided the estimates of r = 0.15 and d r = 60 cm for Ellzey fine sand. Boundary Conditions and Numerical Simulation The governing equations 2 27 and 2 28 are nonlinear partial differential equations whose analytical solutions are not always readily av ailable. For our study, the equations were solved numerically by implicit finite difference method for one dimensional lateral flow. In the field, ditches and water furrows are used to control the WT during irrigation and drainage. The furrows PAGE 42 42 thus act as internal boundaries in additions to the ditches during irrigation and early stage of drainage. Therefore, in the numerical solution, the locations of the water furrows were treated as internal boundaries in addition to the external boundary conditions impo sed at the main ditches. The numerical solution was obtained using the tri diagonal matrix algorithm in the statistical programming language R [ R Development Core Team, 2011]. The simulated WT dynamics were then compared with the field data collected from the study site during the spring of 2010 and 2011.The Nash and Sutcliffe [1970] efficiency coefficient, C eff and root mean square error (RMSE) were used as goodness of fit measures between the model implemented only with d and the model with both d and f Results Drainable Porosity and Fillable Porosity under ET and Recharge Figure 2 7 shows d and f as functions of WT for clay, loam, and sand with typical grass type vegetation at different atmospheric ET demands. Increasing ET demand reduces d to va lues smaller than hydrostatic ( HS ), while at the same WT depths and ET demands, f values are larger than HS This deviation from the HS is different for different soils and is maximum for 3). Fo r example, in loam soil at ET = 0.02 cm/hr, when the WT is near 80cm from the surface, d value is 22% less than the HS while at the same WT depth f is 51% more than HS For clay soil, this difference is not as extreme as in loam while for sandy soils: the difference seems negligible even at ET rate of 0.02cm/hr. As the WT depth increases beyond the ET transition depth, both d and f converge towards HS since the press ure distribution in the soil profile moves towards hydrostatic. When the WT is very close to the surface, both d and f are close to zero irrespective of the ET demand. For sandy soils, which normally have low capillary retention, d and f hardly differ from HS even when the ET demand is high. On the other hand, clay and loam soils show PAGE 43 43 significant differences from HS even under small ET demand because of their high moisture retention capacity. The soil in our study site also did not show as large a shift in d and f from HS as would clay or loam soils under ET (Figure 2 8). This can be attributed to the sandy nature of the profile and shallow root system of potato plants. Nonetheless, there is still a discernible downward shift in d as ET demand increases. Figure 2 8 also suggests that the difference between d and f reaches as much as 18% of the effective porosity ( s r ) at 45cm WT depth for ET = 0.03cm/hr. Hourly ET rates in the study area can reach as high as 0.07 cm/hr during mid day [ FAWN 2011] which would result in significantly lower d or higher f values. The shift in the values of d and f under shallow water table conditions was more pronounced for recharge (Figure 2 9). Under recharge, for the same WT depth, f becomes smaller than HS as R e increases while the opposite is true for d Since R e tends to saturate the soil profile, extra water needs to be drained in order to achieve the same drop in WT that would have occurred under hydrostatic conditions, resulting in larger d On the other hand, f decreases sharply since recharge reduces the storage deficit [ Chapman and Dressler, 1984; Chapman, 1995] consequently reducing the amount of water required to achieve a unit rise in the water table. It should be noted that using dif ferent ET or R e functions may result in slightly different relations between d and f WT for the same soil and vegetation types. Water Table Dynamics The observed and simulated hourly WT dynamics in Well 1 are shown in Figure 2 10 for 50 day periods in sp ring 2010 and 2011. Similar patterns were also observed in the other two wells (data not shown). Model results are compared for the single parameter model (SPM) with HS ( e quation 2 26) and the dual parameter models (DPM) with separate steady state d and f expressions developed in this study ( e quation s 2 27 and 2 28). Note that the ET and R e values PAGE 44 44 required for the estimation of d and f from equations 2 23 and 2 24 were determined a priori which enabled analytical calculation of these storage paramete rs as functions of the predicted WT elevations during each simulation. In both years, the WT data clearly show the timing of irrigation and drainage in the field. During irrigation periods, the WT was consistently maintained around 40 60cm depth from the s urface by raising water levels in the lateral ditches and irrigation furrows accordingly. Diurnal fluctuations in the WT are also clearly discernible during the days with irrigation. This diurnal fluctuation, despite the continuous subsurface lateral inflo w of water, indicates that the WT supports most of the ET demand of potato crop under this water management system. The observed WT dynamics were matched better by the DPM simulations than the SPM. Note that during the early part of the 2011 simulations, t he observed diurnal WT fluctuations were larger than predicted for both SPM and DPM. This indicates that ET for the potato crop during this period was under estimated, likely because the crop coefficient for the mid growth stage of potato in this area may be higher than the general FAO [1998] value. The average DPM C eff values in both 2010 and 2011 (0.56, 0.55) were higher compared to those for SPM (0.34, 0.18), indicating improved performance for the DPM, with an especially marked difference in 2011 (Tabl e 2 2). The average RMSE (cm) values in both years were correspondingly smaller for the DPM (0.34, 0.85) than the SPM (3.14, 0.86). Note that the C eff values for Well 3 were significantly smaller than the values for Well 1 and Well 2. This was likely due t o the proximity of this well to the ditch (12 m) where the assumption of purely horizontal flow may not have been supported due to flow convergence. The SPM simulated WT fluctuated at higher elevations than the observed WT with narrower amplitudes both in 2010 and 2011. The observed diurnal amplitudes of WT depth were well matched by the DPM. This discrepancy between DPM and SPM was more pronounced PAGE 45 45 towards the end of the simulation period both in 2010 and 2011. This period coincided with later stages of th e potato crop when the water level in the ditches and furrows were brought to the shallowest elevation (approximately 30 35cm below ground surface). Rising WT during rainfall events was also predicted by DPM better than SPM, which underestimated WT rise du ring smaller events. These results highlight the contribution of the recharge dependent f parameter in improving the simulation results. Discussions When flux dependent d and f are considered in groundwater studies, it is also important to understand the relationship between these parameters and the magnitude of flux. Of particular interest is the relationship between R e and f because it determines the rate of WT rise during ra infall events. The relationship between f (equation 2 24) and R e (scaled to K s ) under Ellzey fine sand for four different WT depths is shown in Figure 2 11. These results suggest that when the WT is deep, even a small R e corresponds to a f value very clo se to zero. That is, by the time the wetting front has reached a deep WT, most of the soil profile is already saturated. Whereas under the same R e in a shallow WT, some portion of the soil pores may still be empty (i.e. f > 0 ). However, for larger rechar ge rates f approaches zero irrespective of the WT depth. Chapman and Dressler [1984] showed a similar relationship between R e and an effective storage coefficient. An important advantage a method to estimate f in the presence of vertical fluxes is that i t enables simple estimation of WT rise during rainfall events. In shallow phreatic aquifers, the WT normally shows a quick response to the rainfall events. Using f and neglecting the contribution of groundwater flow and ET in equation 2 27, the water tabl e position at time t during rainfall be estimated simply as PAGE 46 46 ( 2 30) and ( 2 31) where t is the time step increment, and P (t ) is the rainfall during the period of each time step. Note that R e is estimated from P using equations 2 19 and 2 20, which is then used to estimate f using equation 2 24. The WT position is then estimated by simple Euler integration of equation 2 30 as s hown in equation 2 31. Application of equation 2 30 to the field site using f is compared in Figure 2 12 to the estimates based on HS The WT rise is significantly underestimated by HS as compared to our f expression, especially when the rainfall rat es are smaller. During high intensity rainfall events, however, the error of estimation is relatively large for both expressions. Comparison of equation 2 31 with HYDRUS 1D for an initial WT depth of 120cm (Figure 2 13) shows good agreement for all rainfal l rates (0.05cm/hour 2.0cm/hr), although equation 2 26 estimates faster WT rise during the initial phase. This is an anticipated pattern because in deeper water tables R e will not occur instantaneously, since the unsaturated zone captures most of the inf iltrated moisture [ Hilberts et al., 2007]. Application of f (equation 2 24) in deeper initial WT conditions therefore may require additional modification of the recharge function to account for the lag between rainfall and recharge. Conclusions Drainable and fillable porosity parameters are critical in modeling groundwater flow in unconfined aquifers. They not only determine the water table dynamics but also influence the subsurface discharge at the aquifer outlet. In this paper we have presented a method to incorporate effects of unsaturated zone moisture fluxes while estimating these two parameters. By assuming successive steady state ET and R e fluxes, closed form expressions for d and f PAGE 47 47 were developed. The main advantage of these expressions is that the dynamic behavior of soil water is accounted for in the estimation of storage parameters. These expressions therefore relax the limitations of previously available methods that relied on a hydrostatic assumption to incorporate the effect of soil moisture retention. The strongest merit of these expressions, however, is the ability to estimate and d and f separately for WT drawdown and rise scenarios. This study shows that depending on the magnitude of vertical flux and soil hydraulic properties, the difference between f and d can be substantial. With the help of a hydraulic groundwater model implemented with flux dependent d and f it is shown that the simulation of WT dynamics can be strongly improved. The values of drainable and fillable p orosity are different under hydrodynamic conditions, but they collapse into a single equation under a static moisture profile. The analytical expressions revealed that under the influence of ET dynamic state d and f of sandy soils hardly differ from HS whereas they can be significantly different in clay and loam. In general, the results suggest at a given WT depth, the finer the soil, the more pronounced the difference between HS and d or f will be for the same ET rate. Recharge also had strong effects on the magnitudes of both d and f It was found that the same recharge rates correspond to a much smaller f values (very close to zero) under deeper WT than under shallower WT. This study h ighlights the importance of using fillable porosity for modeling water table dynamics. The inequality between d and f that develops during recharge or ET in shallow WT environments is significant and should not be neglected in most soil types. Fillable p orosity is normally regarded to be smaller than d due to hysteresis in soil moisture retention [ Bouwer and Jackson, 1974; Bouwer, 1978; Sophoclelous, 1991]. We found differences between d and f even when PAGE 48 48 hysteresis is not considered. Extending these an alyses to incorporate hysteretic effects precludes analytical solutions and requires numerical methods [ e.g., Stauffer et al., 1992]. Note that we examined shallow water table dynamics, but the expressions developed in this study may be effectively used in modeling shallow subsurface flow in hillslopes [e.g., Troch et al., 2003; Hilberts et al., 2004, 2005, 2007], land drainage [e.g. Skaggs 1980], and estimating ET from diurnal water table fluctuations [ Loheide et al., 2005; Nachabe et al., 2005; Gribovszki et al., 2007]. PAGE 49 49 Table 2 1. Original and Modified van Genuchten (VG) soil moisture retention model parameters and s alpha ( a G ) for Clay, Loam and Sand, and Ellzey fine sand series. Original VG parameters or Sand, Clay and Loam are from Carsel and Parris [1988]; a G parameters for Clay, Loam and Sand are from Reynolds and Elrick [1987, 1985]. Parameter Clay Loam Sand Ellzey Original Modified Original Modified Original Modified Original Modified s 0.380 0.398 0.430 0.542 0.430 0.447 0.395 0.398 r 0.060 0.239 0.078 0.087 0.045 0.045 0.066 0.075 0.008 0.000 0.036 0.001 0.145 0.079 0.019 0.011 N 1.09 0.43 1.56 0.47 2.69 1.67 2.63 2.07 a G 0.001 0.0335 0.25 0.068 K s (cm/hr) 0.2 1.04 29.7 7.0 Table 2 2 Nash Sutcliffe Coefficient of Efficiency ( C eff ) and Root Mean Square Error ( RMSE ) of the model implemented with hydrostatic (common d and f ) and dynamic (different d and f ) expressions 2010 Hydrostatic d Dynamic d f C eff RMSE C eff RMSE (cm) Well 1 0.28 3.30 0.80 0.39 Well 2 0.37 3.14 0.75 0.22 Well 3 0.02 2.99 0.12 0.41 2011 C eff RMSE C eff RMSE (cm) Well 1 0.21 0.59 0.75 2.06 Well 2 0.10 1.51 0.63 0.03 Well 3 0.65 2.50 0.28 0.46 PAGE 50 50 Figure 2 1 Conceptual representation of the water table profile and water balance components in a typical unconfined aquifer with shallow phreatic surface, and bounded by a ditch or stream. Figure 2 2. Conceptual representations of change storage deficit ( D s ) due to a unit rise in WT elevation (left), and changes in unsaturated storage ( U s ) and saturated storage ( S s ) due to a unit decline in WT(right) in an unconfined aquifer system PAGE 51 51 Figure 2 3 Relationship between ET rate and WT depth for clay, loam and sand soils under grass type vegetation, based on equations 2 17 and 2 18 The parameter values for the soils were clay: = 95cm, = 0.011; loam: = 85cm, = 0.028; and sand: = 30cm, = 0.043 [Shah et al., 2007]. Figure 2 4. Suction at the soil sur face ( as a function of WT depth calculated from HYDRUS 1D and the analytical expressions ( e quations 2 16 and 2 21) based on the successive steady state assumption during a rainfall rate of 0.03cm/hr for initial WT depth at 50cm (top left) and 75cm (top right), and during WT drawdown due to ET for Ellzey fine sand from the field site. PAGE 52 52 Figure 2 5 Study site in northeast Florida showing agricultural land use and typical water table management system with ditches and irrigation furrows. PAGE 53 53 Figure 2 6 Original and modified van Genuchten moisture retention models fitted to the measured date for Ellzey fine sand series PAGE 54 54 Figure 2 7 Drainable (left) and fillable (right) porosity as the function of WT depth for clay, sand and loam soils under grass type vegetation during hydrostatic and different steady ET rates from WT table based on equations 2 23 and 2 24. The ET parameters were ad opted from Shah et al. [2007]. PAGE 55 55 Figure 2 8 Drainable (left) and fillable (right) porosity of the Ellzey fine sand series soil under potato at different ET demands. PAGE 56 56 Figure 2 9 Drainable (left) and fillable (right) porosity of the Ellzey fine sand series soil under potato for different rainfall (P) rates. PAGE 57 57 Figure 2 10 Simulated and observed hourly water table dynamics in the Well 1 located at the center of the stud y field in 2010 (top) and 2011 (bottom) spring potato seasons PAGE 58 58 Figure 2 11 Plot showing the relationship between recharge rate and fillable porosity for Ellzey fine sand. The relationship is presented in terms log( f ) against to enhance the visual pattern in the curves. PAGE 59 59 Figure 2 12 Observed water table rise at the study site during four rainfall events (3/28, 4/25, 5/04, and 5/17) in 2010, compared to predictions from equations 2 30 and 2 31, neglecting groundwater outflow to the ditches. PAGE 60 60 Figure 2 13 Water table rise predictions during five hypothetical rainfall events in Ellzey soil with initial WT = 120 cm. Lines are from equations 2 30 and 2 31 using recharge dependent f F illed circles are fro m the numerical simulation of the 1D Richards equation using HYDRUS 1D PAGE 61 61 CHAPTER 3 ESTIMATION OF EVAPOTRANSPIRATION FROM DIURNAL WATER TABLE FLUCTUATIONS USING DRAINABLE AND FILLABLE POROSITY Background In shallow groundwater environments such as riparian zone and wetlands, evapotranspiration is a significant component of the water balance [ Dolan et al., 1984; Loheide et al., 2005; Luo and Sophocleous, 2010] where the shallow WT normally supports most of the plant water requirement [ Shah et al., 2007; Grib ovszki et al., 2007; Lowry and Loheide, 2010]. Consequently a distinct diurnal fluctuation pattern in the water table (WT) hydrograph is usually observed [Schilling, 2007; Butler et al., 2007; Nachabe et al., 2005; Loheide, 2008]. This diurnal water table f luctuation (DWTF) corresponds to the 24 hour variation cycles in the climatic variables, solar radiation, temperature and humidity [ Gribovszki et al., 2007; Gribovszki et al., 2010]. A typical DWTF pattern observed in a crop field in northeast Florida is i llustrated in Figure 3 1. Note that the observed WT elevation decreases more rapidly during daytime and slowly resumes during the nighttime. Several different methods are available for estimating the ET from vegetated surfaces which normally involve the measurements of various weather parameters [e.g. Monteith, 1965; Priestley and Taylor, 1972, etc]. The measurements are used to solve different relationships describing the energy balance and aerodynamic vapor transfer phenomena to estimate the evaporative loss. The Penman Monteith method [ Monteith, 1965] is based on the combination of both energy balance and aerodynamic method [ Chow et al., 1964]. One of the main drawbacks of the weather based ET estimation methods is that they usually involve complex calc ulation methods and require large number of parameters. PAGE 62 62 In unconfined aquifers shallow phreatic surface, there is normally a direct link between the vegetation water uptake and the water table (WT). Therefore, ET loss from these shallow WT environments ma y also be estimated from the soil water balance which can be expressed as ( 3 1) where, h is the water table elevation, is the storage parameter of the aquifer, q is the groundwater flux, R e is the rainfall recharge and ET is the evapotranspiration loss. The term in the left hand side of equation 3 1 gives the total change in soil moisture storage per unit area expressed in term of WT elevation. equation 1 can be reformulated to estimate the ETg as ( 3 2) From equation 3 1 it is therefore possible to obtain the ETg if major portion of the ET is directly contributed by the water table. White [1932] was probably the first to analyze and report such direct link between ET and WT fluctuations in shallow WT areas. Based on the observations DWTF in a shallow riparian zone, he used the soil water balance to develop a method to estimate the groundwater water consumption by the phreatophytic vegetation. His method however was based on different sim plifying assumptions. He assumed that rate of groundwater flow to the area was constant over the 24 hour period and could be estimated by the rate of water table rise between midnight and 4:00 am when ETg is negligible or zero. Based on these assumptions a nd ignoring the precipitation component, equation is expressed as ) ( 3 3) where, d is the drainable porosity (or specific yield) of the aquifer, r is the rate of WT rise between midnight and 4:00 AM when the ETg is co nsidered zero, and is the rise or fall of WT during the 24 hour period. In addition to these assumptions, the drainable porosity parameter, d was also assumed as constant in equation 3 3. PAGE 63 63 Following White many have investigated the DWTF phenomenon to estimate the shallow groundwater consumption by vegetation [e.g., Dolan et al., 1984; Gerla, 1992; Loheide et al., 2005, Nachabe et al., 2005, Lautz, 2008, Butler et al., 2007, Gribovszki et al., 2007; Loheide, 2008; Mould et al., 2010, etc.]. Some of these studies also evaluated the underlying assumptions in equation 3 3 and found that the assumptions of constant groundwater flow as well as constant d are normally not met in field conditions [ Gerla 1992 ; Loheide et al., 2005; Gr ibovszki et al., 2007]. Gerla [1992] reported that assuming a constant d in equation resulted in excessively large ETg estimation in a wetland environment. According to Loheide et al. [2005], the assumption of constant d creates problems in the White met hod and showed with the help of soil moisture flow simulations that d was critically important in accurately estimating ETg from DWTF. Loheide [2008] proposed a modified White method to estimate ETg at subdaily time scales considering the variable groundw ater flow rates. Gribovszki et al. [2007] presented a similar modification to Loheide [2008] by using variable, successive steady state groundwater flow and a WT dependent d expression developed by Nachabe [2005]. Although the implementation of variable groundwater flow and variable d have generally shown to resolve the inherent limitations of the White method and allow for subdaily (e.g. hourly) scale ETg estimations [ Gribovszki et al., 2007; Loheide, 2008], there are still some difficulties associated with this method. Firstly, they are based on a d parameter estimated either from the hydrostatic soil moisture profile [ Gribovszki et al., 2007] or from the analysis of rainfall induced WT rise [e.g. Dolan et al., 2007; Schiling, 2007] while this parameter dynamically depends on both WT elevation and soil moisture flux above the phreatic surface [ Childs, 1960; Tritscher et al., 2000]. Since d estimated from hydrostatic moisture profile ( hydrostatic d ) does not account for the un saturated zone flux directly induced by both ET and rainfall infiltration, its PAGE 64 64 value may not accurately represent the actual d in the field under hydrodynamic conditions [e.g. Tritscher et al., 2000]. This discrepancy therefore can introduce potential err ors in the ETg calculations from DWTF. Besides neglecting the effect of soil moisture flux, hydrostatic d expressions also have another drawback; these methods result in the same parameter value for both rising and falling water table conditions while in reality the rise in WT elevation is dictated by a distinct storage parameter called the fillable porosity ( f ). Fillable porosity is an analogous concept to d which is defined as the amount of water absorbed by the aquifer per unit rise in WT elevation [ Bouwer, 1978; Sophocleous, 1991; Park and Parker, 2008]. This parameter is therefore a measure of the volume of pore spaces that remain to be filled during the rise of the phreatic surface [ Chapman and Dressler, 1984; Chapman, 1995]. In Chapter 1 it is sho wn that when vertical flux from and to the water table are considered, d and f tend to differ significantly during ET and recharge thus limiting the use of hydrostatic d in either groundwater flow models or in ETg estimation from DWTF. Since the hydros tatic d ignore s the unsaturated zone soil moisture flux, they cannot accurately reflect rapid WT fluctuation shown by shallow phreatic aquifers during ET and precipitation E specially during and after rainfall events, the assumption of hydrostatic moistur e profile is severely disturbed because of rapid redistribution of soil moisture. Therefore, the actual d may differ significantly from the estimated hydrostatic d [ Tritscher et al., 2000] consequently resulting in potentially erroneous ETg estimations. In order to avoid these complications during rainfall, most of the White based ETg estimation methods presented in the literature typically omitted the days with precipitation in their calculations [e.g. Dolan et al., 1984; Gribovszki et al., 2007; Schilling, 2007]. As a result, these estimations are usually PAGE 65 65 discontinuous at smaller time scales (i.e. hourly or daily) although fairly well estimates may be obtained on larger time scales (e.g. monthly average) if the rainfall frequency is low. W hile the potential impact of rainfall on the underlying assumptions of White based ETg estimation methods results due to the disturbance in the DWTF signal of the WT, the literature provides very little work on the extent to which the ETg estimates are aff ected when days with rainfall are included in the calculation. Since the method is based on soil water balance, it can be argued that this method should work satisfactorily even during the rainfall events, if the effect of rainfall on WT fluctuations can b e quantified reliably. In one study, Gerla [1992] included the days with precipitation to estimate ET using the White method. He modified the original White method by accounting for the effect of infiltration on water table rise thus extending the use of t he method even during precipitation events. His study however lacked an appropriate analytical expression to estimate the storage parameters and used d values estimated from the WT rise during infiltration. In addition, method was based on daily t ime scale where the hourly dynamics of recharge and WT rise are lumped thus concealing the potential impacts of rainfall infiltration on the hourly ETg estimates. Nonetheless, Gerla [1992] showed a significant improvement over other the original White meth od, since it could be used during the rainfall events to obtain the ETg estimates that correlated well with the potential ET (PET) estimates from weather measurements. The main purposes of this paper is therefore twofold: (1) To use a d parameter estimat ed from the soil moisture profiles under successive steady state soil moisture flux ( dynamic d ), instead of the hydrostatic d to estimate ETg from DWTF and compare the results from the two methods, (2) To include the days with rainfall also in the calcu lations by modifying the White method to account for the effect of precipitation on the WT fluctuation using the f parameter PAGE 66 66 also estimated under successive steady soil moisture flux. The dynamic d and f partially incorporate the effect of the soil moi sture flux on the respective changes in aquifer storage and deficit during WT movement due to ET and/or groundwater flux. They may provide potentially more reliable approximation of the actual storage parameters of the aquifer than given by the hydrostatic d while also enabling the application of appropriate parameter s to WT ri se and drawdown respectively. Theoretically, use of these distinct d and f parameters should enable the estimation of ETg from soil moisture balance continuously regardless of the rainfall events because it is possible to estimate the contribution of rainfall to change in WT. Therefore, an attempt is made in this study to obtain the continuous estimates of ETg without avoidi ng rainy days in order to comparatively assess the impact of rainfall on ETg estimates. Firstly, the governing equations to estimate ETg implemented with corresponding d or f parameter are presented based on whether the WT is being driven downward or upw ards due to the groundwater flow. The method of estimating the two storage parameters are also briefly discussed. A simple method of attributing the WT rise to recharge during rainfall based on the WT hydrograph and f is also illustrated. This step provid es a scheme to extract the exact WT change due to the groundwater flow and ETg if there was no precipitation during a particular time step. Finally, the equations are applied to estimate ETg from a potato field under ditch drainage and irrigation system in northeast Florida. Both hydrostatic d parameter method and the new method are used to estimate hourly and daily ETg over a 50 day period during spring 2010 and 2011 seasons. The results from the two methods are then compared with each other with respect to the standard Penman Monteith method [ Monteith, 1965) as a benchmark, estimated from the measurements of weather parameters in the study site. PAGE 67 67 Theory Estimating ETg from DWTF using d and f The soil water balance given by equation 3 2 provides the total change in storage expressed in terms of WT elevation measured from a reference datum. The parameter in equation 3 2 can be either d or f depending on whether the WT is receding or rising as suggested by their respective definition s [ Bouwer, 1978; Sophocleous, 1991]. In order to implement these two parameters in the water balance equation it is therefore necessary to understand the respective effect of each water balance component on the water table. For example, ET drives the water table downwards from the soil surface. So, the total change in aquifer storage caused by ET can be estimated by using d Recharge ( R e ) on the other hand, normally causes the phreatic surface to rise, which is dictated by the storage deficit of the aquif er and hence by f The groundwater flow ( q ) can result in either rise or drawdown WT depending on the boundary conditions thus requiring appropriate storage parameter to determine the total change in storage. Therefore, based on whether the water table is rising or falling due to groundwater flow, the water balance equation can be reformulated appropriately by implementing d and f as ( 3 4) during groundwater inflow with recharge ET and ( 3 5) during groundwater outflow with recharge and ET Note that in equations 3 5, the R e term is associated with both d and f This is because in a declining WT scenario, the total amount of recharge at the phreatic surface is partitioned into two fractions. A fraction that contributes to recharge and cause the water table to rise, and PAGE 68 68 another fraction that is drained with the g roundwater outflow and lost from the profile. Therefore the effect on WT is determined by the difference between the reciprocal of f and d Equations 3 4 and 3 5 can now be reorganized to give the governing expressions implemented with distinct d and f to estimate ETg as ( 3 6) for ETg from WT during groundwater inflow (rising WT) and recharge. ( 3 7) for ETg from WT during groundwater outflow (draining WT) and recharge In equation 3 6 and 3 7, the term represents the contribution of recharge towards the rise in WT elevation, so it needs to be accounted for and removed from the observed change in WT elevation if a particular day with recharge is to be included in ETg calculations.This implies that the resulting equation is not applicable to days with recharge when the water table is declining. Implementing separate d and f parameters thus allows to include days with rainfall with improved ETg estimations and generate continuous subdaily ETg calculations. In addition, immediately after the rainfall event which causes the water table rise, ET loss normally occurs at the atmospheric potential rate. This results in d values that are smaller than the hydrostatic d consequently estimating the WT downward more quickly than the latter This eventually produces smaller ET estimates from equation 3 7 than from equation 3 2. Estimation of d and f Application of equations 3 6 and 3 7 to estimate ETg necessi t ates a method to estimate d and f as the function of WT elevation and unsatured zone moisture fluxes which was presented in Chapter 2. These d and f expressions partially account for these unsaturated zone water flow PAGE 69 69 by assuming steady state flow at successive time steps. The generalized expressions for d and f under hydrodynamic conditions in the soil profile given are given as ( 3 8) ( 3 9) where, s = saturated water content; r = residual water content; = matric suction at the soil surface; and and n are the parameters of a modified van Genuchten [1980] soil moisture retention model [ Troch, 1992; Bierkens, 1998; Hilberts et al., 2005 ]. The value of can be estimated by using hydaulic conductivity function [ Gardner, 1958] assuming steady state flows. Note that, both equation 3 8 and 3 9 collapse into the same expression when hydrostatic conditions are assumed in the soil profile (i.e., The magnitudes of d and f given by equation 3 8 and 3 9 therefore depend on the soil hydraulic parameters, the magnitudes of ET or R e and the vegetaion types. Figure 3 2 shows the drainable and fillable porosity of Ellzey sand series soil from northeast Florida under potato duri n g different rates of ET and recharge respectively. The plots suggest that the difference between d and f incre ases as the magnitude of the fluxes increase. Note d and f for ET are plotted under relatively large hourly ET values becase it has been found that the ET loss in the area may reach as high as 0.07cm/hr during the midday [ FAWN, 2011). Athough the differ ences in d and f for the Ellzey soil during ET are not as high as they woul be in finer textureed soils, it is still substantial to be replaced by a single, hydrostatic parameter. Hourly Water Table Slopes Normally during rainless days, hourly changes in the WT elevation at each time step is simply estimated by substracting the observed WT elevation at a previous time step from the PAGE 70 70 observation at the current time step. However, whenever there is rain fall, becomes large due to WT recharge, which eventually results in highly negative ETg estimations. It is due to this complication that most of the previous studies on ETg have typically avoided the days with substantial precipitation [e.g. Dolan et al., 1984; Gribovszki et al., 2007; Schilling, 2007]. However, if estimates of R e and f are available, the contribution of recharge on the observed time series of can be approximately quantified and removed from the WT time series [ Healy and Cook, 2002). Removal of the recharge induced WT rise at each time step then yields the estimate of the actual change in WT that would have been observed if there was no rainfall during that period. This approach is fundamentally simila r to the methods that estimate recharge from WT time series [e.g., Crosbie et al., 2005] by extracting the WT level rise that can be attributed solely to recharge [ Healy and Cook, 2002; Crosbie et al., 2005]. However, in our study, instead of estimating r echarge from the time series as done by Crosbie et al. [2005), we use a simple exponential function to estimate R e based on the WT elevation and rainfall, which upon division by the f yields time series. Figure 3 3 shows the schematic illustration of the sudden rise of an initially declining WT during a single time step (dt) whose magnitude is given by line BD. The observed value of is the the ratio AB/BD which reflects the net effect of both recharge and the groundwater outflow Note that in this particular case the groundwater flow drives the water table downward while the recharge drives it upward. Line DC in Figure 3 3 represents the anticipated locus of the WT elevation if there was no recharge [ Healy and Cook, 2002). This a nticapated path of WT can be determined by estimating the groundwater flux at point D and dividing this value by d The ratio PAGE 71 71 AC/BD is the contribution of recharge on the WT rise. The exact change in WT which is the ratio BC/BD, is then simply obtained as the difference between AC/BD and AB/BD. Note that when the water table recession is very slow during the time period from point D to C, i.e., if shallow u nconfined aquifers of large lateral later extent where the groundwater flow is much slower as compared to the vertical recharge during rainfall. Extraction of from requires the a reliable estimation f parameter [ Healy and Cook, 2002] which changes dynamically during the infiltration process [ Chapman and Dressler, 1984; Chapman, 1995]. The method employed in this study to extract can be summarized as follows Firstly, the observed values are calculated from the entire dataset which is the total effect of ETg q and R e The next step requires estimation of R e from any appropriate method which then enables the calculation of f as the function of the WT elevation and R e During recharge, the amount of water that accretes on the WT may be much larger than the amount removed/added via groundwater flow. Therefore, neglecting the groundwater flow during a particular time step with recharge, can be estimated as ( 3 10) Note that in equation 3 10, in addition to neglecting the groundwater flow, instantaneous recharge is assumed, which is reasonable for shallow water tables [ Hilberts et al., 2007; Laio et al., 2009]. However, this assumptions may not hold in deeper WT environments w here a lag in WT response is likely to occur during precipitation events as the infiltrated water need longer time to reach the water table. This lag in WT response, if likely to occur, can be determined by comparing the time series of rainfall and the values. If a lag in relative to the time of PAGE 72 72 rainfall is observed, i.e, if perceptible response of WT to R e is not detected in that particular time step, then the calculated from equation 3 10 is applied to the following time step that shows the effect of recharge (by sudden rise in WT). After this process is repeated over the entire time seris of WT elevation, can be simply calculated as ( 3 11) During intensive rainfall events, estimated values become very high, even yielding infinite values since large recharge rates drive f towards zero [ Chapman and Dressler, 1984]. In such cases, can be replaced by the observed s ince will not exceed the observed by a large magnitude at a given time step, unless the groundwater flux rate is very high. Figure 3 4 shows the observed and the actual calculated during three separate rainfall events at a field site in northeast Florida. ETg Calculations Once the effect of precipation have been removed from the observed time series, equations 3 6 and 3 7 can be solved to obtain the ETg estimates as follows: during groundwater inflow, ( 3 12) ( 3 13) and during groundwater o utflow, ( 3 14) ( 3 15) PAGE 73 73 where, t denotes the time step (hourly in this study) and all other symbols are as previously defined. Note that in equation 3 12 the term is set equal to unity. This was because the original equation produced highly unstable ETg estimates as compared to equations 3 12 and 3 13. This was possibly because of the dynamic nature of equations 3 12 and 3 13 since they require estimations of d and f using the ETg value calculated at the previous time step to compute ETg on the current time step. Moreover, if the magnitude of q in the aquifer is small, the first term on the right in equation 3 12 and 3 13 will not change by much even when the is not treated as unity. Also, in equation 3 14 the recharge term is replaced by the product of f and which gives the actual value of recharge that was estimated during the calculation of if the method involved accounting for the lag in WT response. This procedure therefore will ensure consistency in the overall calculation procedures. Application of the Modified Method Site Description and Water Table Data The field site chosen for the application of modified White method developed in this study is located in northeast Florida. This area is typically characterized by low, flat landscape with with a shallow water table predominantly planted to potato. The fields in the area are irrigated by a conventional water table In this system, shallow water table is maintained near the crop root zone during the entire cropping season [ Campbell et al., 1978] by simultaneously supplying water from individual water furrows and raisin g the water level in the ditches Water is then supplied to the crop root zone by capillarity and subsurface lateral flow [ Pitts and Smajstrla, 1989; Smajstrla et al., 2000]. Figure 3 5 shows location of the study sites near the Atlantic coastline in north east Florida. Each PAGE 74 74 seepage irrigated field is divided into approximately 18 20m wide crop beds separated by smaller, shallow drains (or furrows) approximately 0.30 0.50m wide. During irrigation the drains are continuously supplied with water which brings W T near the root zone. The level of water in the drains and ditches thus determines the ultimate height of WT in the field. It can be hypothesized that almost all of the ET demand of the crops is contributed by the shallow WT under this water management sys tem. It is based on the results of several previous studies on shallow groundwater areas [e.g., Nachabe et al., 2005; Shah et al., 2007; Lowry and Loheide, 2010 etc.]. In our study site, WT is very close to the surface which implies a more plausible scenar io of full contribution of WT to plant ET. These field sites therefore also provide a setting to test this hypothesis by estimating ETg from observed diurnal water table fluctuations. The soil at the sites is classified as Ellzey fine sand series (sandy, s iliceous, hyperthermic, arenic endoaqualf; 90 95% sand, <2.5 % clay, <5% silt) [ NRCS, 1999]. Although this series is predominated by sand, Acharya and Mylavarapu [2011] studied different physical properties of the soil at the study site at varying depths a nd reported that the clay content in the profile increased from 2% at 0 22cm to 11% at 90 120 cm. They also reported that the lab determined average saturated hydraulic conductivity of the soil profile from 0 120 cm depth was approximately 10 cm / hr while in an in situ measurement performed by Rosa [2000], it was found that K s of the soil was approximately 7 cm/hour. The latter was used in this study because in situ measurement is more likely to represent the real world situation as compare d to the lab dete rmined values. Two monitoring wells, 1.1 m deep were installed in two crop beds near the center of a potato field of approximately 15ha (400m x 400m) area in the spring of 2010. Hourly water table elevation in the wells was recorded using the pressure tran sducers for 50 days from March to PAGE 75 75 May. In spring 2011, in addition to the field 1, two wells were installed in a second potato field (field 2) of similar dimensions and the water table elevations were recorded for approximately 50 days. Water table elevat ions recorded by the pressure transducers were then adjusted for the variations in barometric pressure recorded at the center wells of each field to obtain the corrected WT elevations. The approximation positions of the the wells in the two field sites is shown in Figure 3 5. Estimation of Recharge The recharge values were estimated using an empirical R e Vs.WT function which depends both on rainfall rate and WT elevation and is determined by the soil type. As shown in Chapter 2, for estimating the rise in shallow WT during rainfall events, this functions approximated the WT rise resonably well in the stud y site. The empirical function, which depends on the soil type can be expressed as ( 3 16) ( 3 17) where, P is the rainfall amount during the given time step, and r and d r are the parameters that determine what fraction of rainfall is converted to recharge depending on the soil type. Equation 3 16 and 3 17 thus provide a method of estimating recharge and the c onsequent WT rise during rainfall without explicitly considering the unsaturated zone transport processes [ Woods et al., 1997]. The values of d r and r which were determined by fitting the equations with the WT rise obtained from hypothetical infiltration simulations by Hydrus 1D [ Simunek et al., 1998], were estimated as 60 cm and 0.15 respectively. All other parameters required and used in this study were adopted from the earlier study conducted in the field sites. PAGE 76 76 Estimation of Subsuface Flux ( q ) In area s with WT control systems, ditches and drains are appropriately constructed to raise or lower the WT elevations as needed during the crop growht period. Subsurface flow in these fields is dictated by the water level in the drains which can be treated as en tirely horizontal, adopting the Dupuit Forcheimmer (D F) theory [ van Schilfgaarde, 1974; Ritzema, 1994]. Figure 3 6 shows a conceptual representation of the WT profile in the study field with main lateral ditches and smaller inner ditches (or, furrows) dur ing an irrigation period. During the drainage phase the water table profiles are exactly opposite (i.e. convex WT profile) than during irrigation phase. From the D F theory, if an instantaneous steady state flow at each time step is assumed, the subsurfac e lateral flux to and from the ditches which fully penetrate to the impervious layer can be estimated as (Ritzema, 1994) ( 3 18) where, = flux towards the outer ditches; K s = lateral saturated hydraulic conductivity of the soil; = height of WT at midpoint between the ditches, measured from impermeable layer; = height of the water level in the ditches; and L D is the distance between the main ditches. However, if the d rainage ditches are shallower, the flow converges near the drains which needs to be accounted for while estimating the subsurface flux. Flow to and from the shallow and ( 3 19) where, m = difference between the height of WT at the drains and WT at midpoint between the drains; L d = the distance between the smaller drains which is normally 18m in the study fields. and d e = equivalent depth of th e impermeable layer introduced to account for the convergence of PAGE 77 77 flow near the shallow drains [ Moody, 1966; van Schilfgaarde, 1974; Ritzema, 1994]. Moody [1966] derived a simple formula to calculate the equivalent depth of the impermeable layer from drain radius, drain, depth of the impervious layer, and distance between the drains, which was used in this study. However, any other appropriate method [e.g., van Beers, 1979; van der Molen and Wesseling, 1991] may be used to estimate the equivalnet depth. Onc e both and are estimated, the total flux can be simply obtained as Penman Monteith ET Estimation of ET from Penman Monteith (PM) method [ Monteith, 1965] requires the data on solar radiation, temperature, relative humidity, and wind speed. These weather data were obtained from a Florida Automated Weather Network (FAWN) weather station adjacent to the field sites [ FAWN, 2011]. The standardize PM meth od for hourly reference ET [ Allen et al., 1998] is expressed as ( 3 20) where ET 0 = reference evapotranspiration (mm/hr); R n = net solar radiation (MJ/m 2 h); G is the soil heat flux (MJ/m 2 hr); U 2 = wind = slope of the saturation vapor pressure curve (kPa/C) at mean air temperature (T); e s = saturation vapor pressure at the mean hourly temperature (kPa); e a = acutal vapor pressure at mean hourly temperature; = latant heat of vaporization (MJ/kg); C d = bulk surface and aerodynamic resistance coefficient. Equation 3 20 was used to calculate the hourly reference ET ( ET 0 ) for the field site during the spring o f 2010 and 2011. The reference evapotranspiration was then converted to the actual ET for potato using the potato crop coefficient reported by Allen et al. [1998]. The estimations were obtained for approximately 50 days in both years starting from the mid growth stage of potato until one or two weeks before harvesting. PAGE 78 78 Results and Discussions Hourly ETg Estimates Hourly estimates of ETg from the governing equation implemented only with hydrostatic ( HSM, hereafter ) and dynamic d and f ( DM hereafter ) compared with the standardized PM estimates during the entire study periods in 2010 and 2011 are presented in Figure 3 7. The data presented in Figure 3 7 are from the well located at the center of field 1 in 2010 and field 2 in 2011. However, the data f rom the other wells also showed highly similar results in both fields (data not shown). Note that during both seasons the days with rainfall were included in the calculations of ETg by both HSM and DM Generally during rainless days, the HSM significantly overestimated ET as compared to the DM Overestimation of ETg by the HSM was more prominient in field 1 during 2010 than in field 2 during 2011. The generally higher amplitudes of the DWTF in 2010 than in 2011 (Figure 3 7) might have caused s uch higher degree of overestimation in 2010. Note that in 2010, during majority of the study period, average amplitude of diurnal WT fluctuations was more than 10cm while in 2011 it was less than 10cm. This was because of the generally low potential ET during the 2011 study period than during 2010. The weather data from the weather station revealed that, although the daily average temperatures during 2010 and 2011 in the site did not differ significantly (data not shown), 2010 was relatively drier than 2 011 during the study period. This generally higher relative humidity migh have played a role in producing higher PET values during 2010 which, on average, was 0.1cm more than in 2011 study period (March May). The average RMSE of hourly ETg estimates in 201 0 and 2011 seasons were 0.012 cm and and 0.0008cm for the HSM and DM respectively suggesting that the DM was able to vastly improve the ETg estimates from the field sites during the study periods. PAGE 79 79 Since one of our main objectives was to study how the Whi te method based hourly and daily ETg estimates are affected during and immediately after precipitation events, the results from HSM and DM models were compared specifically during rainy and clear days. Figures 3 8 and 3 9 show the hourly ETg estimates from HSM DM and PM ET superimposed on the associated WT elevation curve during and immediately after two rainfall events and two rainless periods in 2010 and 2011. Note that even on the days with precipitation, hourly pattern of ETg estimates did not deviat ed greatly from the normal diurnal pattern (i.e. as shown by PM ET). It was found that during the rainfall events also, the ETg estimates from DM were fairly close to the PM ET while HSM tended to overestimate. In general, there was a significant difference in the ETg estimates form HSM and DM immediately after rainfall when the WT was rapidly receding due to drainage and ET. However, for the rainfall event on 05/05 in 2010 and during the WT drawdown period follo wing the event (Figure 3 8), the two method were mostly similar to each other and to the PM ET. The PM ET, on the other hand, showed relatively smaller values during this 3 day period as compared to the days before rainfall (05/04) and the times after WT r ecession stopped (05/09). This was also reflected by the relatively less steep WT recession curve during that period as compared to the other drawdown events (e.g. 4/26 4/27 in 2010 and 5/16 5/18 in 2011). These observations are consistent with the theory of flux dependent d which suggests that as the flux rate decrease, the difference between hydrostatic d and dynamic d decreases and vice versa. As ET rates become smaller, the corresponding soil moisture profile above WT moves towards the hydrostatic profile and the two parameters become essentially the same As a result DM also behaves as the HSM producing similar ETg estimates. These results suggest that the method used in our study to PAGE 80 80 account for contribution of recharge in observed perfor med reasonably well, improving the ETg estimates during and after most rainfall events throughout the two study periods. A closer evaluation of hourly ETg estimates during the WT drawdown periods after each rainfall show that the result from DM model tend to fluctuate substantially near the end of each drawdown event, i.e., when the WT recession pattern halts as the irrigation is resumed. For example in Figures 3 8 and 3 9, a significant fluctuations in the daytime ETg estimates from DM is observed on 4/15 and 4/28 in 2010 and 4/10 and 4/19 in 2011. This was due to the associated lag in the WT response to the irrigation from the smaller ditches (furrows). In the field, as the requires a few hours in order for the water in the furrows to build up,infiltrate towards the WT, and begin lateral movement. However, in our calculations, this lag in WT response to irrigation was not considered since the infiltration from the smaller di tches (furrow) was not modeled. Therefore, during these particular times, there is a discripancy between the observed WT slope and the associated groudnwater flux in the calculation s, resulting in the fluctuating ETg values. Daily ETg Estimates The daily v alues of ET and estimated ETg also followed similar patterns observed in the hourly estimates. Averag e daily PM ET for the potatoes in 2010 and 2011 were approximartely 0.65 and 0.55cm while the average daily ETg estimates were 0.95 & 0.82 for HSM and 0.64 and 0.56 for DM respectively. Daily ETg estimates from HSM methods were also significantly higher than both DM and PM estimates during the clear days which is presented in Figure 3 10. Average RMSE of daily estimates were 0.28 cm and 0.009 cm for the twe methods in 2010 and 2011 respectively, indicating a massive improvement in ETg estimation by DM over the HSM method. PAGE 81 81 The potential PM ET values estimated in our study fall within th e range of the previously reported potential ET values for the region [e.g., Jones et al., 1984; Smajstrla et al., 1985]. However, data on continuous crop coefficient ( K c ) of potato in the area is still lacking which implies some degree of uncertainty on t he actual PM ET values estimated for the potato crop. Majority of our study period included the mid growth stage of potatoes when the water demand is at the highest level. During the mid season growth, the FAO reported K c for potato [ Allen et al., 1998] i s 1.10, which was used in our study. Other studies have suggested that this value may be higher than 1.1 for the area. Singleton [1990] reported monthly K c values as high as 1.35 for the region based on the time of planting which would have produced actual PM ET values in the range of 0.65 0.75 for our study period. However, due to the lack of continuous data on K c we used the FAO values which are generally regarded as the standard crop coefficients. Neverthless, it should be noted that even if the actual potato K c values were higher, the average daily ETg estimates from HSM method would still be much higher than the PM ET as compared to the DM suggesting the better performance by DM Overall, the DM significantly improved the ETg estimates by HSM method d uring majority of the study periods in 2010 as well as 2011. Interestingly, the daily ETg estimates on the rainy days were essentially similar for both HSM and DM at almost all instances (Figure 3 10). One of the primary reasons behind this may be the oc currence of mostly nighttime precipitation. During nighttime highly negative ETg estimates were converted to zero since during this time period ET may be negligibly small in most cases. Another cause may be that even the daytime rainfalls generally lasted only for a few hours in most precipitation events. As shown in Figures 3 7, 3 8 and 3 9, the longest rainfall durations were only a few hours during bot h 2010 and 2011 study periods. After each rainfall, the WT recession curves show a clear effect of diur nal ET pattern. That means the disturbances PAGE 82 82 to diurnal WT patterns due to rainfall lasted only briefly. Since the HSM method generally overestimated ETg during immediate hours after the rainfall events, the 24 hour sum during those days was still similar to the DM estimates as well as the PM ET. On the contrary to the exact days with rainfall, the HSM method significantly overestimated both DM an PM ET wh ich was not surprising since the hourly estimates showed clear overestimation on most of the rainless days. In both 2010 and 2011, there were only a few days on which the DM method overesti m ated PM ET. The daily sums therefore also corroborate the hourly r esults suggesting that our method performed satisfactorily to handle the rainy days during ETg estimations from the diurnal WT fluctuations. The most critical aspect of estimating ETg from DWTF, attributing the observed appropriately to ETg which is controlled by the d Our study showed that during clear days, the ETg dependent drainable porosi t yt ( dynamic d ) predicted hourly estimates significantly better than the hydrostatic d which consistently produced higher values. Because the unsaturated zone soil moisture flux has a substantial effect in the magnitude of drainable porosity [ Childs, 1960; Tritscher et al., 2000], the WT drawdown is more appropriately attributed to ETg when dynamic d This effect is not incorporated in hydrostatic d so it attains a greater value than the dynamic d for the same WT elevation resulting in larger ETg values. This is also explained by the observed dynamics of the two parameters during the ETg computation process. Figure 3 11 shows the hourly variations in hydrostatic d and dynamic d due to the WT fluctuation during WT recession after rainfall and during clear days in 2010. Note that at daytime, dynamic d values decrease by as much as 50% of hydrostatic d values while at the nighttime, ET falls close t o zero and dynamic d collapses into the hydrostatic d PAGE 83 83 An important aspect of the method presented in this study is that ETg and the storage parameters are estimated dynamically by using estimated ETg values from the the previous time step. This process sometimes may cause some fluctuations in the d and f and hence in the hourly ETg estimates. However, the overall results were generally consistent througho ut the study periods in both 2010 and 2011 which strongly indicate that the method can be used ef fectively to estimate ETg form the diurnal water table fluctuations. Conclusions In this study we have presented an application of distinct, flux dependent d and f parameters to estimate subdaily as well as daily ETg usin g a modified form of the White method. Most of the previous studied based on this method used either constant d [ Lautz, 2008; Mould et al., 2010] or a single d estimated either from static soil moisture profile [ Gribovszki et al., 2007] or from analysis of WT rise during rainfall eve nts [ Dolan et al., 1984; Schilling, 2007]. Use of these separate, flux dependent storage parameters enable the modification of the White method by accounting for the effect of rainfall on the observed WT fluctuations. This modification therefore allows fo r the reasonable application of the White method even during rainy days. In order to include the rainy days in the calculations, a simple method to extract the recharge induced WT fluctuation from the observed hourly changes in WT was also employed in our study. A critical aspect of this method is the appropriate estimation of WT recharge from rainfall which also affects the magnitude of f The empirical recharge function used in the study produced satisfactory results for the shallow WT in northeast Flori da potato fields providing good estimates of the recharge induced WT rise. For deeper WT, however, the method may need slight modifications to account for the lag in the WT response. Our method showed good estimations of the Penman Monteith ET for potato in during spring of 2010 and 2011 even when precipitation events were included in the calculations. PAGE 84 84 Drainable porosity plays a critical role in estimating ETg from White based method s since this method attempts to attribute the observed WT fluctuations to direct ET loss from groundwater via this d parameter. Our study also showed that the d estimated from hydrostatic soil moisture profile normally failed to appropriately attribute the WT fluctuations to the groundwater ET, producing significant overestim ations. However, when the effect of ET is incorporated into the d its value become smaller than hydrostatic d since this dynamic moisture profile based parameter also accounts for the soil moisture that is being lost from the profile via ET, in addit ion to the water lost by drainage. These smaller d values therefore resulted in more reliable estimates of ETg at hourly as well as daily time scales. Our study was based on the assumptions of non hysteretic soil moisture retention properties. It is well known that hysteresis also tends to cause significantly different d and f values [ Nachabe, 2002; Nachabe et al., 2004; Loheide et al., 2005]. Since these hysteretic parameters are applicable in falling and rising water tables respectively, the resulting equation of ETg will also need modification if hysteresis is considered [e.g. Loheide et al., 2005]. However, in field conditions the h ysteretic effects may not always be significant due to other sources of variations. The results from our study further accentuate the critical role of d and f in White based ETg estimation methods while providing a general framework of applying the meth od continuously regardless of the precipitation events using these two storage parameters. Although our conclusions are based on relatively short term results from a managed crop field, close agreement between the results from two consecutive years suggest that the method is consistent and can be equally applied to other field conditions. PAGE 85 85 Table 3 1 Root Mean Square Error (RMSE) of hourly and daily ETg estimations using single, hydrostatic d and distict, dynamic d and f RMSE (cm) Hourly Daily 2010 Eq d 0.013 0.29 NonEq d and f 0.0007 0.009 2011 Eq d 0.012 0.27 NonEq d and f 0.001 0.009 PAGE 86 86 Figure 3 1. Diurnal fluctuations in the solar radiation (red) and observed water table depth (blue) at a northeast Florida site Figure 3 2. Drainable and fillable porosity of Ellzey sand series under potato in northeast Florida under different potential ET and precipitation rates PAGE 87 87 Figure 3 3. Illustration of contribution of recharge induced sudden WT rise during a rainfall even t and extraction of the exact change from the observed by subtracting the contribution of recharge Figure 3 4. Calculated hourly actual change in water table elevation from the observed data after removing the contr ibution of recharge during four rainfall events in spring 2010 PAGE 88 88 Figure 3 5. Land use of the three counties in northeast Florida and the two field sites with well positions, chosen for the study Figure 3 6. Con c eptual representation of the watertable profiles and the flow towards the ditches in the field PAGE 89 89 Figure 3 7. Time series of the hourly rainfall, Penman Monteith and groundwater ET predicted by using a single drainable porosity (middle plots) and distinct drainable and fillable porosity (bottom plots) from field site 1 in 2010 (top), and field site 2 in 2011 (bottom) during the spring potato seasons PAGE 90 90 Figure 3 8. Hourly ETg from the HSM method (black dashed line), DM method (red solid line), and PM ET (blue dots) during clear periods (top) and two rainfall events in 2010. Figure 3 9. Same as Figure 8, during 2011 in field site 2. PAGE 91 91 Figure 3 10. Scatter plots of daily Penman Monteith ET and predicted ETg from field site 1 in 2010 (left) and field site 2 in 2011 (right) during spring potato seasons Figure 3 11 Plots showing the variations in hydrostatic d (open circles) and dynamic d (filled diamonds) during WT dra wdown after rainfall (top) and diru nal fluctuations in the WT elevation in a field site in 2010 PAGE 92 92 CHAPTER 4 MODELING SHALLOW WATER TABLE DYNAMICS UNDER A SUBSURFACE IRRIGATION AND DRAINAGE SYSTEM IN FLORIDA Background Shallow water table (WT) commonly exists in many areas of the world where field hydrology of such differs significantly from well drained soils with deep water table. The hydrology as well as growth and development of plants in such SWT areas are controlled by the dynamics of the phreatic surface [ Bierkens, 1998; Nac habe et al., 2002]. Agricultural crop production in areas with shallow WT therefore requires an efficient drainage system to create optimum environment from crop growth and field operations [ Schilfgaarde, 1974; Tang and Skaggs, 1977]. Drainage is used to r emove excess soil moisture from the root zone [ Fipps and Skaggs, 1989]; to control potential salinity problems in irrigated arid areas [ Fipps and Skaggs, 1989; Ritzema 1994]; and to discreetly manage the water resources. In humid areas such as Florida, sin ce the groundwater level is close to the surface, even a small amount of precipitation may bring WT close to the surface [ Gillham, 1984; Novakowski and Gillham, 1988; Hilberts et al., 2007]. This requires rapid drainage to avoid excessive moisture stress t o plants Such rapid WT drawdown is normally enabled by constructing ditches and installing drains in the field. However, since the drains and ditches hasten the drainage process, the soil around the crop root zone is likely to dry quickly. Therefore many o f these subsurface drainage systems are operated as WT management systems which can be used for controlled, subsurface drainage as well as irrigation whenever necessary In Florida, groundwater level is normally close to the surface in majority of the agr icultural areas C rop production in these areas has been conventionally performed by managing the WT throughout the growing seasons. Water table is maintained at a desired elevation near the crop root zone using open ditches and appropriately spaced surfac e water PAGE 93 93 furrows, conventionally known as seepage irrigation system. Moisture supply to the plants in this system occurs by subsurface lateral flow and capillarity [ Pitts and Smajstrla, 1989; Smajstrla et al., 2000]. Since the WT is very close to the root z one, majority of the crop evapotranspiration (ET) demand is contributed by the direct upward flow from the water table [ Nachabe et al., 2005]. The knowledge of WT dynamics due to diurnal ET forcing, frequent rainfall events, and water level changes in fie ld ditches is therefore important for maintaining optimum soil moisture around the crop root zone. The process of WT movement due to subsurface drainage/irrigation or rainfall/ET from the soil profile is governed by the principles of water movement through soils. Therefore, the design and analysis of such systems can typically be done with the help of phenomenon in soils [ Pikul et al., 1974; Tang and Skaggs, 1977; Hilbe rts et al., 2005]. However due to the inherent difficulties in numerical solutions of this equation, approximate methods that describe the drainage process are normally adopted [ Tang and Skaggs, 1977]. Most of the research on the subsurface drainage of agr icultural lands and hillslopes are based on the Dupuit Forcheimmer (D F) assumptions to describe groundwater flow in unconfined aquifers [e.g, van Schilfgaarde, 1963; Moody, 1966: Skaggs, 1991; Verhoest and Troch, 2000; Rupp and Selker, 2005; Singh and Jaiswal, 2006; etc]. The D F assumption states that in shallow, unconfined aquifers with small WT inclinations, the vertical hydraulic gradient can be neglected and the flow can be assumed to occur entirely in the horizontal direction. Adoption of these as sumptions results in Boussinesq type groundwater flow models, also known as the hydraulic groundwater model [ Brutsaert, 2005]. PAGE 94 94 The Boussinesq type groundwater flow models consist of a storage parameter which enables the expression of aquifer storage in te rms of water table elevation. T he response of WT to groundwater inflow and outflow is therefore controlled by this parameter. When the water table is receding, the parameter is called as the drainable porosity ( d ) which gives the change in aquifer storage per unit decline in WT depth [ Bouwer, 1978; Freeze and Cherry, 1979]. In case of WT rise the parameter is termed as fillable porosity ( f ) which is controlled by the moisture deficit of the aquifer [ Bouwer, 1978; Healy and Cook, 2002; Park and Parker, 20 08]. It is now well known that the moisture storage, and hence WT dynamics of unconfined aquifers is significantly affected by the capillary properties of the soils [ Gillham, 1984; Nielsen and Perrochet, 2000; Healy and Cook, 2002, etc.]. Some authors have incorporated the effect of soil moisture retention in estimation of aquifer storage parameters to study subsurface flow and water table dynamics in hillslopes and wetland environments [ Hilberts et al., 2005, 2007; Laio et al ., 2009]. However, all of the previous studies on the water table response to subsurface drainage and irrigation in managed agricultural lands treated this parameter as constant [ Skaggs, 1991, Singh and Jaiswal, 2006 etc.]. This can introduce substantial e rrors in estimated WT rise or drawdown rates since the effect of unsaturated flow due to ET and recharge is not incorporated in parameter estimation. As shown in Chapter 1, although treated as a single common parameter estimated from static soil moisture p rofiles in recent studies [ Bierkens, 1998; Hilberts et al., 2005], these two parameters are affected by the unsaturated fluxes above the phreatic surface and can assume substantially different values depending on the magnitude of ET or recharge fluxes. In this study we numerically solved the 1D Boussinesq equation implemented with separate, flux dependent d and f parameters in contrast to the conventional Boussinesq model PAGE 95 95 with only one, hydrostatic d parameter The dual parameter model is therefore theoretically more robust since it appropriately implements the storage parameters depending on the effect of various water budget components on the direction of WT movement. The main objective of the study was to eval uate the effectiveness of the WT management system to maintain the desired WT elevation throughout the cropping period to create favorable soil moisture environment around the crop root zone. The second objective was to evaluate the response of WT to poten tial alternative management scenarios. Firstly conceptual models of WT management system under subsurface irrigation and drainage system are presented. These WT management systems may consist of a typical ditch drainage design as well as more complex syste ms with smaller, surface or subsurface drains in between the two main ditches. Appropriate boundary conditions are then applied to the system during irrigations and drainage to solve the governing equations by implicit finite difference method. Simulated w ater table elevations are then compared with the field data collected from a potato fields in northeast Florida managed under a subsurface irrigation/drainage system. The sensitivity of the water table movement to changes in soil parameters such as K s and the field characteristics, and potential alternative water management system specific to northeast Florida region is briefly discussed. Materials and Methods Study Site The field site chosen for this study is located in northeast Florida ( near the Atlantic coastline (Figure 4 1) The flat and low (<10 masl) landscape with with a shallow groundwater level in the area therefore requires appropriate water table management measures. This is achieved by constructing larger (1 2 m deep) ditc hes around the field as well as smaller and shallow surface drains conventionally know n as Each seepage irrigated field is divided into approximately 18 20 m wide crop beds separated by smaller drains PAGE 96 96 (or furrows) approximately 0.30 0 .50 m wide. The shallow water table is maintained near the crop root zone during the entire cropping season [ Campbell et al., 1978] by supplying water from the inner drains and raising the water level in the ditches simultaneously Moisture is supplied to the crop root zone by capillary action and subsurface lateral flow [ Pitts and Smajstrla, 1989; Smajstrla et al., 2000, Munoz et al., 2008; Acharya and Mylavarapu, 2011]. The water supply in the smaller drains is cut off during and after rainfall vents and the water level in the ditches is lowered thus allowing for quicker water table drawdown. The ultimate water table elevation in the crop beds is thus determined by the water level in the ditches and the drains which change during irrigation and drainage. T he soil at the field sites is classified as Ellzey fine sand series (sandy, siliceous, hyperthermic, arenic endoaqualf; 90 95% sand, <2.5 % clay, <5% silt) [ NRCS, 1999]. Although this series is predominated by fine sand, soil contains some clay content which ranged from 2% at shallower depth to 11% at 120 cm [Acharya and Mylavarapu, 2011). Lab determined averave vertical saturated hydraulic conductiivty (K) of the soil was 10 cm/hr while the insitu determined K values of the soil was found to be 7cm/hou r [Rosa, 2000]. Model Development Conceptual models: irrigation and drainage Figures 4 2 and 4 3 show the conceptual representation of the water table profile during irrigation and drainage respectively in the study fields. At the start of an irrigation ev ent, water is supplied to the inner drains while the water level is raised in the ditches from an initial height ( h 0 ). Note that in Figure 4 2, initial water table profile is horizontal while in drained lands the water table normally tends to occur as curv ed profile [ Skaggs, 1973; van Schilfgaarde, 1974]. We assumed a horizontal profile because of the large distance between the two ditches (normally PAGE 97 97 > 250 m) due to which the water table profile is usually close to horizontal near the end of a drainage event Whenever there is an anticipated rainfall event, the water supply to the field is cut off, and water level in the ditches is lowered thus triggering subsurface drainage of the profile. Normally after rainfall, the drainage occurs in two distinct phases as shown in Figure 4 3. If the water table is raised close to the surface (i.e., above the depth of the furrows), both the outer ditches as well as the furrows contribute towards the drainage. This helps is rapid drawdown of the water table during the init ial hours after cessation of rainfall. This is termed as the phase 1 drainage in Figure 4 3. As the water table moves below the level of the depth of the furrows, they cannot contribute to drainage anymore and the system is converted to a typical ditch dra inage system with only the outer ditch removing the subsurface water from the profile. This is termed as phase two drainage in Figure 4 3. Phase 2 drainage is normally much slower than the phase 1 drainage due to large distance between the two ditches. Gov erning equations Since the subsurface drainage phenomena occur largely by the lateral movement of soil moisture below the water table, it can be described by the one dimensional Boussinesq equation subject to recharge and evapotranspiration. Specific gover ning equations implemented with the drainable and fillable porosity for the dynamics of the water table are discussed in Chapter 1. The equations for the drainage and irrigation phases differ only on the use of the storage parameter and can be expressed as ( 4 1) during drainage when the water table is receding ( 4 2) PAGE 98 98 during irrigation when the water table is rising, and where, h = the ele vation of water table measured from a reference datum K s = saturated hydraulic conductivity of the profile d = drainable porosity f = fillable porosity ET = evapotranspiration R e = Recharge Note that in equation 4 1 and 4 2, the both d and f appear instead of one common parameter as used in the previous studies. The processes that are governed by these equations are discussed in detail in Chapter 1. In order to solve 4 1, and 4 2, d and f should be estimated as the function of h by incorpor ating the effect of unsaturated zone flux since ET and R can largely affect their magnitudes [ Childs, 1960, Chapman and Dressler, 1984, Chapman, 1995; Tritscher et al., 2000]. As derived in Chapter 1, specific expressions for f and d during ET and R e are given as ( 4 3) ( 4 4) Where = the matric suction at the soil surface = saturated water content profile = residual soil water content of the profile and n = parameters of the modified van Genuchten [1980] soil moisture retention equation [ Troch, 1992; Hilberts et al., 2005] PAGE 99 99 In equations 4 3 and 4 4, the soil matric suction at the surface is estimated from the water Gardner, 1958] and assuming successive steady state vertical flow to and from the water table [ Raats, 1974; Zhu and Mohanty, 2004] as shown in Chapter 1. Initial and boundary conditions Irrigation : The initial and boundary condition for the irrigation are similar to those for the phase 1 drainage (equation 4 5). Once the water supply to the furrows is initiated, water infiltrates vertically towards the WT Since the width of the furrows is very small as compared to the field width, they can be considered as line sources [ Gill, 1984]. After the wetting front PAGE 100 100 re aches the WT, the soil profile below the furrow s is completely saturated and therefore can be considered as narrow vertical ditches from which water seeps laterally. If the WT is deep (e.g., > 2m), it may require a long time for the wetting from to reach t he water table and contribute to WT rise. In our study fields, however, the wetting front quickly reaches the phreatic surface due to high saturated hydraulic conductivity of the soil and shallow WT depth. Therefore, it was assumed during the solution that the soil profile below the furrow is saturated immediatel y after the onset of irrigation. N umerical solution Equations 4 1 and 4 2 are nonlinear partial differential equation which cannot be solved analytically for many field geometries and numerical meth ods are generally adopted. In this study we solved the governing equations numerically by implicit finite difference method for the initial and boundary conditions expressed in equations 4 5 and 4 6. Equation 4 1 can be expressed in finite difference form as ( 4 7) S etting equation 4 7 can be, rearranged to get a system of linear equations with three unknowns, expressed as ( 4 8) where, ; ; ; & Eq uation 4 8 was solved by using the Tri diagonal matrix algorithm using R programming environment [ R Development Core Team, 2012] to estimate hourly water table profiles within the field. PAGE 101 101 Water Table Data and Soil Parameters Water table monitoring wells eac h approximately 1.2m deep, were installed at different locations in the field site in 2010. Figure 4 1 shows the location of the wells in the field which are located between the mid point of the field and one of the ditches. Hourly water table data was rec orded from the wells during the spring potato season of 2010 using pressure transducers. In the following year, hourly water table elevations were recorded in field site 1 as well as field site 2, which is approximately of same dimensions (400 m x 400 m). The data recorded by the pressure transducers were compensated to remove the effect of barometric pressure in the observed water table elevations. Water levels in the water furrows and ditches, which provide the boundary conditions for the water flow in th e field, were treated as constants during the irrigation and drainage events thus creating constant head type boundary conditions. The height of the water in the furrows during irrigation and drainage phases normally differed by5 10cm. The beginnings of ir rigation and drainage phases were obtained by recording the time of turning on/off of the water supply in the field. Hourly rainfalls as well as other weather parameters (solar radiation, temperature, relative humidity, wind speed, and dew point temperatu re) were obtained from the weather station install adjacent to the study fields [ FAWN, 2011]. Evapotranspiration from the potato fields were calculated at hourly time steps from the Penman Monteith method [ Monteith, 1965] using the hourly weather data. To calculate the recharge rate from rainfall data, a parametric relationship was used which is discussed in Chapter 1. All the parameter required for the simulations, along with the procedures of their estimation, are presented in Chapter 1. Simulation of wa ter table dynamics during irrigation and drainage events were performed on hourly time scales for 50 days period during the mid growth season of potatoes in 2010. The simulated water table elevations at the monitoring wells were then compared to the field water PAGE 102 102 table data. In 2011 potato season, the model was tested against the two independent datasets from field 1 as well as field 2. The two statistical measures, Root Mean Square Error (RMSE) and Nash Sutcliffe Coefficient of efficiency ( C eff ) [ Nash and Su tcliffe, 1972] were used to assess the performance of the simulation. Once the performance of the WT dynamics simulation model was evaluated from the two year potato season data, the sensitivity of the model to few key parameters and alternative irrigation regimes was studied. One of the major drawbacks of the seepage irrigation system is that it is highly inefficient in terms of water use since it requires continuous pumping of water to the furrows during irrigation. Therefore, a water saving irrigation re gime, with only 12 hour water supply (7:00 am 7:00 pm) was formulated the resultant water table dynamics was studied with respect to the dynamics observed under conventional irrigation regime. Results and Discussions Model Application The statistical measu res used to compare the simulated and observed water table elevations during the 50 day simulation periods in 2010 and 2011 are presented in Table 4 1. Figure 4 4 shows the water table depths against time (hours) in Well 1 during two drainage/rainfall even ts and two irrigation events in 2010. Note that there is a significant difference in the water table movement rates during irrigation and drainage phases which are due to the effect of ET which affects drainable and fillable porosity resulting in a hystere tic behavior of the two storage parameters (Chapter 2), which eventually causes significant discrepancy in water table movement. During the irrigation phase, even though the water level in the furrows (internal boundary condition) was normally maintained a round 40 cm, the WT fluctuated diurnally at 50 60cm depth showing the effect of direct ET from WT The depth of the root zone of potato crop in the area is reported to be approximately 36cm [ Munoz et al., 2006]. These PAGE 103 103 observations therefore indicate that the conventional system is effective in consistently maintaining WT close the bottom of the root zone thus avoiding any moisture stress during irrigation. Simulated rates of WT rise due to rainfall as well as the subseq uent recession rates matched closely with the observed data in all three locations. The best agreement between the observed and simulated water table dynamics over the 50 day period was found in Well 1 which was located at the center of the field. The aver age C eff and RMSE (cm) of the simulations in 2010 were 0.56 and 0.51 respectively with Well 3 showing the lowest C eff This suggests that the model was not able to capture the water table dynamics as well in the periphery of the field as it did near t he ce nter. One of the primary reasons for this is t he underlying assumption of the model itself where exclusively horizontal flow is assumed following the D F theory. N ear the ditches, the streamlines tend to converge instead remaining parallel [ van Schilfgaar de, 1974; Ritzema, 1994] which can overestimate the WT movement. Figure 4 5 shows the scatter plots of observed of simulated WT depth in the three wells in 2010 which corroborate the results of Table 4 1. The disagreement in observed and simulated WT depth occurred mainly during heavy rainfall events. Th e model overestimated the WT drawdown especially after larger rainfall events potentially because of the effect of surface storage. The surface storage can infiltrate towards the WT after the rainfall ceases thus slowing down the drawdown process. Such slower drawdown rates can also be detected in Figure 4 4 during the two precipitation events. Model Evaluation and Sensitivity Analysis The performance of the WT dynamics model was also tested against two indep endent WT datasets collected from field 1 and fiel 2 during the 2011 potato season. Figure 4 6 shows the plots of observed versus simulated water table depths in field 2 during the 50 day period in 2011 The statistical measures of comparison between the s imulated and observed data for both fields PAGE 104 104 are presented in Table 4 1, which strongly support the results from 2010 potato season. The average C eff and RMSE (cm) of simulations in 2011 were 0.55 and 0.85 for field 1, and 0.71 and 0.56 respectively. These comparisons suggest that the numerical model developed in this study reliably described the water flow processes that govern the drainage and irrigation in the field sites. Although the numerical Boussinesq model presented in our study was able to simulate hourly water table dynamics during both irrigation as well as drainage phases in the field sites, it is critical to understand the sensitivity of the model outcomes to various parameters. One of the most critical parameters that determine the rate of WT r ise and drawdown in the field is the saturated hydraulic conductivity ( K s ) Figure 4 7 presents the effect of different K s rates on WT drawdown and rise after the onset of the drainage and irrigation respectively Water table movement highly depend s on the K s of the soil since this parameter determines the rate of soil water flow. In our particular fields, when the K s value increased from 2 to 7 cm/hr, the time required to bring the WT initially at 20cm depth, below 40cm decreased dramatically to < 24 hour s form approximately 50 hours. Similarly, when the K s value increased to 15 cm/hr the drawdown time further decreased approximately to 12 hours. Water table rise during irrigation event, like the drawdown phase, also differed significantly when the K s valu es changed. The time required for the WT, initially at 80cm below the surface to reach near the bottom of the root zone (40cm) was < 24 hours for K s = 15cm/hr while it was 100 hours for K s = 7. 0 cm/hr. The WT hydrograph observed in the field after the ons et of the irrigation events also showed the same pattern takin g approximately 100 hours to reach the maximum WT elevation where it fluctuated diurnally due to the influence of ET. The maximum WT elevation reached during the irrigation event was also higher with greater K s values. These results can have an important PAGE 105 105 implication in irrigation scheduling under conventional WT management system in Florida since the irrigation onset times can be optimized based on the anticipated WT movement rates I n addition to the K s of the soil, a critical component of subsurface drainage system is the depth and spacing of drains. In most parts of the Florida including our field sites, narrow water furrows are created approximately at 18 m distance between the mai n lateral ditches which facilitate both irrigation and drainage. Therefore, it is important to understand their contribution to water table drawdown. Figure 4 8 presents the effect of furrow depth on the water table recession below the potato crop root zon e in a field that resembles the field sites used in this study. Note that increasing the drain depth just by 10 cm can significantly reduce the time required to bring WT below the root zone. When the furrow depth is increased to 65cm from 45cm, which is th e normal depth of the furrows practiced in northeast Florida, time required to bring water table below 40cm was reduced by more than 20 hours. Also note that after the water table passes the depth of the furrow, the drawdown curve shows a nearly linear rel ationship to time suggesting that the contribution of ditches to water table drawdown is small compared to the furrows in these fields. The recession of WT after it passes the depth of the furrows therefore mostly might have occurred due to the direct ET f rom the phreatic surface [e.g., Loheide et al., 2005; Shah et al., 2007] which attenuates the effect of the ditch water level on subsurface lateral flow towards the field edges. A cri tical limitation of continuously maintaining the WT near the root zone is that it requires continuous pumping of water, a substantial portion of which is lost to the ditches via runoff from furrow. Due to diminishing groundwater resources, it is necessary to optimize the water management system to reduce water loss and increase the efficiency. A potential PAGE 106 106 alternative to the conventional system would be an intermittent system in which water is supplied only during the day time instead of the night time. It i s important to know how the water table responds to the intermittent water supply system in order to understand the potential impact on crop growth. Figure 4 9 shows the simulated water table dynamics under a hypothetical intermittent irrigation system as compared to the observed dynamics during two irrigation events in 2010 potato season. Note that under the intermittent system, the diurnal pattern of the water table is still similar to the continuous system while the maximum elevation of the water table i s decreased by as high as >15cm. This suggests that the intermittent system cannot maintain the water table as close to the root zone as the continuous system does, which may potentially induce moisture stress to the root near the soil surface. However, fo r crops with relatively deeper root system, this system may be effective enough. Since the intermittent system save 50% of water application, its potential as an alternative water management system can be justified. However, extensive research on the root zone soil moisture dynamics is necessary to understand the feasibility of such water management regimes in the area. Conclusions In this study we used numerical technique to solve the one dimensional Boussinesq equation to model the shallow water table dynamics in unconfined aquifers and applied to in a subsurface irrigation and drainage system typically practiced in Florida. Our mo del was developed by implementing two distinct drainable and fillable porosity parameters which respectively control the water table drawdown and rise in unconfined aquifers. This approach is significantly different from all the previous studies in which o nly a single drainable porosity was used. Our model was able to simulate the recession and rise of water table during respective drainage and irrigation events in northeast Florida. Direct ET loss from the phreatic surface, which commonly occurs in shallo w water table environment, was captured very well by the numerical model. PAGE 107 107 Application of the model to the field sites in northeast Florida also suggested that, in addition to the soil hydraulic properties, the water table movement during irrigation and dra inage is mostly dictated by the depth of the water furrows and ET while the effect of ditches at the field boundaries was comparatively negligible. This was further supported by the simulated water table responses at varying depths of the water furrows. Si mulated water table dynamics during a water saving, potential alternative to the conventional water table management revealed that the water table remained at approximately 15cm lower depth as compared to the current system. Although this alternative saves 50% of water pumping to the field, the possibility of lower water table elevations may induce moisture stress to the plant roots at shallow depths. Nonetheless, for relatively deeper rooted crops, the results indicate the feasibility of implementing such intermittent system in the area. Subsurface irrigation and drainage system are necessary in shallow water table environments to facilitate crop growth and field traffic while they can also reduce the pollutant loads to surface waters. Numerical models not only help understand various soil water phenomena in the field but they also help in testing various alternative management scenarios that may be potentially adopted to improve the efficacy of current management systems. Since majority of the agricultural lands in Florida are under subsurface irrigation and drainage system, the numerical model presented in this study can be helpful in improving efficiency of the water management system as well as reduce intensive groundwater withdrawal in the region through optimization of irrigation scheduling. PAGE 108 108 Table 4 1 Comparison measures: Root Mean Square Error (RMSE ) and Nash Sutcliffe Coefficient of Efficiency of water table dynamics simulation for the three wells in site 1 in 2010 and 201 1. Distance From the Main Ditch C eff RMSE (cm) Site 1 2010, Field 1 Well 1 200 0.80 0.39 Well 2 125 0.75 0.22 Well 3 13 0.12 0.41 2011, Field 1 Well 1 200 0.75 2.06 Well 2 125 0.63 0.03 Well 3 13 0.28 0.46 2011, Field 2 Well 1 200 0.78 0.54 Well 2 125 0.76 0.12 Well 3 13 0.59 1.04 PAGE 109 109 Figure 4 1. Land use in the St. Johns, Putnam and Flagler counties in northeast Florida and a field site, planted to potato, chosen for the study PAGE 110 110 Figure 4 2. Conceptual representation of water table profiles during irrigation phase. Note that only 2 water furrows are depicted here; the actual number of furrows in the fields may be10 30 Figure 4 3. Conceptual representation of water table profile during the early drainage phase (top) and late drainage phase (bottom). The dashed line represents the initial water table elevation after a precipitation event PAGE 111 111 Figure 4 4. Observed (filled circles) and simulated (solid lines) water table rise due to rainfall (vertical bars) and subsequent drawdown (top), and water table rise after onset of irrigation (bottom) in Well 1 located approximately at the center of the field site 1 in 2010 potato season. PAGE 112 112 Figure 4 5. Scatterplots of the observed and simulated hour ly water table depths in the three wells in 2010. The red line indicated 1:1 relationship. Figure 4 6. Scatterplots of observed and simulated hourly water table depths in field site 2 in 2011 potato season. The red line indicate s 1:1 relationship. PAGE 113 113 F igure 4 7. Effect of K s on the time required to lower the water table below 50 cm depth from an initial WT depth of 10 cm (left) ; and time required to raise the WT to 50 cm from an initial depth of 80 cm (right) in a field of 15 ha area (400m x 400 m). Figure 4 8. Effect of furrow depth s on the water table drawdown rate during subsurface drainage PAGE 114 114 Figure 4 9. Hourly water table depths during irrigation events of an alternative, 50% water saving irrigation regime (12 our irrigation/drainage cycle) PAGE 115 115 CHAPTER 5 SUMMARY AND CONCLUSIONS In C hapter 2, a method to incorporate effects of unsaturated zone moisture fluxes while estimating d and f was developed. By assuming successive steady state ET and R e fluxes, closed form expressions for d and f were derived. The major advantage of the method is that the dynamic behavior of soil water is accounted for in the estimation of storage parameters thus relaxing the limitations of previously available methods that relied on a hydrostatic assumption to i ncorporate the effect of soil moisture retention. Most important advancement, however, is the ability to estimate and d and f separately for WT drawdown and rise scenarios. With the help of a hydraulic groundwater model implemented with flux dependent d and f it is shown that the simulation of WT dynamics can be strongly improved. Our study therefore has highlighted the importance of using fillable porosity which can be critical in modeling WT dynamics in land drainage systems [e.g. Skaggs 1980], hillslopes [e.g., Troch et al., 2003; Hilberts et al., 2004, 2005, 2007] and wetlands [e.g., Laio et al., 2009; Pumo et al., 2010]. In Chapter 3, an application of the two analytical expressions d and f developed in C hapter 2 to ETg from observed diurnal water table fluctuations observed is presented. In doing so, the original method developed by White [1932] and later modified by others [e.g. Gribovszki et al., 2007], was further modified to account for the effect of precipitation in observed WT fluctuations. This modification therefore allowed for the reasonable application of the White method even during rainy days. Our study also showed that the d estimated from hydrostatic soil moisture profile normally failed to appropriately attribute the WT fluctuations to the groundwater ET, producing significant overestimations. However, when the effect of ET is incorporated into the d its value becomes smaller than hydrostatic d since this non equilibrium parameter also accounts for the soil moistu re that is being lost from the profile via PAGE 116 116 ET, in addition to the water lost by drainage. These smaller d values therefore resulted in more reliable estimates of ETg at hourly as well as daily time scales. The results from our study further accentuate the critical role of d and f in White based ETg estimation methods while providing a general framework of applying the method continuously regardless of the precipitation events using these two parameters. Although our conclusions were based on relatively short term results from a managed crop field, c lose agreement between the results from two consecutive years suggest that the method is consistent and can be equally applied to other field conditions. In C hapter 4, a numerical model of WT fluctuation under a subsurface drainage and irrigation system is developed by numerically solving the 1D Boussinesq equations. The governing equations resulted by implementing separate d and f as shown in C hapter 2 were applied in the study. The model was then applied to a subsurface irrigation and drainage system ty pically practiced in northeast Florida. The model simulated the drawdown and rise of water table during respective drainage and irrigation events in northeast Florida. Direct ET loss from the phreatic surface, which commonly occurs in shallow water table e nvironment, was captured very well by the numerical model. Application of the model to the field sites in northeast Florida also suggested that, in addition to the soil hydraulic properties, the water table movement during irrigation and drainage is mostly dictated by the depth of the water furrows and ET while the effect of ditches at the field boundaries was comparatively negligible. Subsurface irrigation and drainage system are necessary in shallow WT environments like northeast Florida to facilitate cro p growth and field traffic while they can also reduce the pollutant loads to surface waters. Numerical models not only help understand various soil water phenomena in the field but they also help in testing various alternative management scenarios that may be potentially adopted to improve the efficacy of current management systems. Therefore, the model presented in this PAGE 117 117 study can be potentially helpful in improving efficiency of the water management system as well as reduce intensive groundwater withdrawal in the region through optimization of irrigation scheduling. The dynamics of shallow WT in unconfined aquifers largely control the hydrology and vegetation growth development in natural areas while it presents a great challenge to efficient water management in agricultural lands. Such dynamic behavior is in turn con trolled by the physical characteristics of the aquifer material as well as the vegetation. Drainable and fillable porosity are the two primary aquifer parameters that control the WT fluctuation in shallow phreatic aquifers. They are critical in modeling gr oundwater flow in unconfined aquifers and not only determine the water table dynamics but also influence the subsurface discharge at the aquifer outlet. This research primarily focused on investigating the influence of unsaturated zone flux above the WT on the behavior of d and f which was not done until now. The results revealed that their relationship with the WT is manifested differently under hydrodynamic conditions. These results thus paved a way to formulate more robust expressions of the hydraulic groundwater model that can be applied in agricultural WT management systems, shallow hillslope subsurface flow, or WT fluctuations in wetlands and riparian zones. Furthermore, the results enabled the investigation of close interactions between shallow WT a nd plant ET through the development of more detailed and reliable calculation methods. Collectively, these improvements aid to better understanding of the complex dynamics of soil plant water interactions in shallow WT environments which are present all ar ound the globe and help solve more problems in the future. PAGE 118 118 APPENDIX THE WATER TABLE DYNAMICS SIMULATION MODEL: R CODE seepageModel < function( pars ) { N = constants $ Time RL = constants $ Z.eq delt = constants $ delt delz = constants $ delz z = seq ( 0 pars $ zmax delz ) L = length ( z ) Bed.Len = pars $ zmax #### Dupuit Forcheimmer Flux Function flux.DF < function( wt ditch.wl ) { return ( pars $ Ksat wt ( wt ditch.wl ) / delz ) } ### Determine number of Furrow/Drains Bed.num < function( pars constants ) { Rd1.wd = constants $ Rd1.wd Rd2.wd = constants $ Rd2.wd Bed.wd = pars $ Bed.wd Fur.wd = pars $ Fur.wd Fur.no < ( pars $ zmax ( Rd1.wd + Rd2.wd ) + pars $ Bed.wd ) / ( pars $ Bed.wd + pars $ Fur.wd ) Bed.no < as.integer ( Fur.no 1 ) return ( Bed.no ) } #### Determine Furrow locations FurLoc < function( pars constants ) { Fur.no = Bed.num ( pars constants ) + 1 xx < rep ( pars $ Bed.wd / delz Fur.no ) xy < 1 : Fur.no 1 Ist.fur = constants $ Rd1.wd / delz + 1 xxz < xx xy xz < Ist.fur + xxz yy < rep ( pars $ Fur.wd / delz length ( xz )) fur < recur ( xz yy ) Furrow.Loc = fur + xy return ( Furrow.Loc ) } ##### Drainable Porosity under ET sdPor.ET < function( pars constants Ht ET ) { e = min ( ET ET exp ( pars $ b ( pars $ Td Ht ))) qv = e psiT < function( Ht qv ) { return ( 1 / pars $ a log ((( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv ) / pars $ Ksat )) } A < function( Ht qv ) { return (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) / (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv )) } PAGE 119 119 sdP < function( psiTop gama pars ) { return ( pmax (( pars $ WCs pars $ WCr ) ( 1 gama ( 1 + ( pars $ alpha siTop ) ^ pars $ n ) ^ ( ( 1 + 1 / pars $ n ))), 0 )) } dpor = sdP ( psiT ( Ht qv ), A ( Ht qv ), pars ) return ( dpor ) } dPor.ET < cmpfun ( sdPor.ET ) ##### Fillable Porosity Under ET sfPor.ET < function( pars constants Ht ET ) { e = min ( ET ET exp ( pars $ b ( pars $ Td Ht ))) qv = e psiT < function( Ht qv ) { return ( 1 / pars $ a log ((( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv ) / pars $ Ksat )) } A < function( Ht qv ) { return (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) / (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv )) } sfPor < function( psiTop gama pars ) { fill.por = ( pars $ WCs pars $ WCr ) gama ( 1 ( 1 + ( pars $ alpha psiTop ) ^ pars $ n ) ^ ( ( 1 + 1 / pars $ n ))) return ( fill.por ) } fpor = sfPor ( psiT ( Ht qv ), A ( Ht qv ), pars ) return ( fpor ) } fPor.ET < cmpfun ( sfPor.ET ) #### Drainable Porosity under Rainfall/Recharge sdPor.Rf < function( pars constants Ht Rf ) { qv = min ( Rf Rf exp ( pars $ br ( pars $ Tr Ht ))) psiT < function( Ht qv ) { return ( 1 / pars $ a log ((( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv ) / pars $ Ksat )) } A < function( Ht qv ) { return (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) / (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv )) } sdP < function( psiTop gama pars ) { return ( pmax (( pars $ WCs pars $ WCr ) ( 1 gama ( 1 + ( pars $ alpha psiTop ) ^ pars $ n ) ^ ( ( 1 + 1 / pars $ n ))), 0 )) PAGE 120 120 } dpor = sdP ( psiT ( Ht qv ), A ( Ht qv ), pars ) return ( dpor ) } dPor.Rf < cmpfun ( sdPor.Rf ) #### Fillable Porosity under Rainfall/Recharge sfPor.Rf < function( pars constants Ht Rf ) { qv = min ( Rf Rf exp ( pars $ br ( pars $ Tr Ht ))) psiT < function( Ht qv ) { return ( 1 / pars $ a log ((( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv ) / pars $ Ksat )) } A < function( Ht qv ) { return (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) / (( qv + pars $ Ksat ) exp ( pars $ a ( Ht RL )) qv )) } sfpor = function( psiTop gama pars ) { fill.por = ( pars $ WCs pars $ WCr ) ( gama ) ( 1 ( 1 + ( pars $ alpha psiTop ) ^ pars $ n ) ^ ( ( 1 + 1 / pars $ n ))) return ( fill.por ) } fpor = sfpor ( psiT ( Ht qv ), A ( Ht qv ), pars ) return ( fpor ) } fPor.Rf < cmpfun ( sfPor.Rf ) #### Thomas Algorithm For Tridiagonal Matrix Tridia < function( A B C R ) { A = as.vector ( A ) ##First Diag B = as.vector ( B ) ## True Diagonal C = as.vector ( C ) ## 3rd Diag R = as.vector ( R ) nn = length ( A ) BETA < rep ( 0 nn ) GAMA < rep ( 0 nn ) BETA [ 1 ] = C [ 1 ]/ B [ 1 ] GAMA [ 1 ] = R [ 1 ]/ B [ 1 ] for ( i in 2 : nn ) { BETA [ i ] = C [ i ]/ ( B [ i ] BETA [ i 1 ] A [ i ] ) GAMA [ i ] = ( R [ i ] GAMA [ i 1 ] A [ i ] ) / ( B [ i ] BETA [ i 1 ] A [ i ] ) } Y < as.vector ( rep ( 0 nn )) Y [ nn ] = GAMA [ nn ] for ( j in 1 : ( nn 1 )) { Y [ nn j ] = GAMA [ nn j ] BETA [ nn j ] Y [ nn j + 1 ] } PAGE 121 121 return ( Y ) } Tridiag < cmpfun ( Tridia ) #### Coefficient Vector Functions: Second Phase Drainage (Ditch Only Scenario) AA.Drain < function( Ht ) { A < as.vector ( rep ( 0 length ( Ht ))) nn = length ( A ) for ( i in 2 : ( nn 1 )) { A [ i ] = Ht [ i 1 ] A [ 1 ] = 0 A [ nn ] = 0 } return ( A ) } AA.D < cmpfun ( AA.Drain ) #### BB.Drain < function( lamda Ht ) { B < as.vector ( rep ( 0 length ( Ht ))) nn = length ( B ) for ( i in 2 : ( nn 1 )) { B [ i ] = 1 / lamda [ i ] + Ht [ i 1 ] + Ht [ i ] } B [ 1 ] = 1 B [ nn ] = 1 return ( B ) } BB.D < cmpfun ( BB.Drain ) #### CC.Drain < function( Ht ) { C < as.vector ( rep ( 0 length ( Ht ))) nn = length ( C ) for ( i in 2 : ( nn 1 )) { C [ i ] = Ht [ i ] } C [ 1 ] = 0 C [ nn ] = 0 return ( C ) } CC.D < cmpfun ( CC.Drain ) #### RR.Drain < function( pars BC1.D BC2.D Ht dp fp lamda ET Rf ) { e = ET Re = min ( Rf Rf exp ( pars $ br ( pars $ Tr Ht ))) R < as.vector ( rep ( 0 length ( Ht ))) nn = length ( R ) for ( i in 1 : nn 1 ) { R [ i ] = ( Ht [ i ] delt e / dp [ i ] + delt Re ( 1 / fp [ i ] 1 / dp [ i ] )) / lamda [ i ] } R [ 1 ] = BC1.D R [ nn ] = BC2.D PAGE 122 122 return ( R ) } RR.D < cmpfun ( RR.Drain ) #### Coefficient Vectors: First Phase Drainage: Ditches and Furrow/drains AA.Drain2 < function( Ht ) { fur.loc < FurLoc ( pars constants ) A < as.vector ( rep ( 0 length ( Ht ))) nn = length ( A ) for ( i in 2 : ( nn 1 )) { A [ i ] = Ht [ i 1 ] A [ 1 ] = 0 A [ fur.loc ] = 0 A [ nn ] = 0 } return ( A ) } AA.D2 < cmpfun ( AA.Drain2 ) #### BB.Drain2 < function( lamda Ht ) { fur.loc < FurLoc ( pars constants ) B < as.vector ( rep ( 0 length ( Ht ))) nn = length ( B ) for ( i in 2 : ( nn 1 )) { B [ i ] = 1 / lamda [ i ] + Ht [ i 1 ] + Ht [ i ] } B [ fur.loc ] = 1 B [ 1 ] = 1 B [ nn ] = 1 return ( B ) } BB.D2 < cmpfun ( BB.Drain2 ) #### CC.Drain2 < function( Ht ) { fur.loc < FurLoc ( pars constants ) C < as.vector ( rep ( 0 length ( Ht ))) nn = length ( C ) for ( i in 2 : ( nn 1 )) { C [ i ] = Ht [ i ] } C [ fur.loc ] = 0 C [ 1 ] = 0 C [ nn ] = 0 return ( C ) } CC.D2 < cmpfun ( CC.Drain2 ) #### RR.Drain2 < function( pars BC1.D BC2.D BC3.D Ht dp fp lamda ET Rf ) { e = ET Re = min ( Rf Rf exp ( pars $ br ( pars $ Tr Ht ))) fur.loc < FurLoc ( pars constants ) R < as.vector ( rep ( 0 length ( Ht ))) PAGE 123 123 nn = length ( R ) for ( i in 2 : nn 1 ) { R [ i ] = ( Ht [ i ] delt e / dp [ i ] + delt Re ( 1 / fp [ i ] 1 / dp [ i ] )) / lamda [ i ] } R [ 1 ] = BC1.D R [ nn ] = BC2.D R [ fur.loc ] = BC3.D return ( R ) } RR.D2 < cmpfun ( RR.Drain2 ) #### Coefficient Vector Functions: Irrigation Case: Same as Ist Phase Drainage AA.Irrig < function( Ht ) { fur.loc < FurLoc ( pars constants ) A < as.vector ( rep ( 0 length ( Ht ))) nn = length ( A ) for ( i in 2 : ( nn 1 )) { A [ i ] = Ht [ i 1 ] A [ 1 ] = 0 A [ fur.loc ] = 0 A [ nn ] = 0 } return ( A ) } AA.I < cmpfun ( AA.Irrig ) #### BB.Irrig < function( lamda Ht ) { fur.loc < FurLoc ( pars constants ) B < as.vector ( rep ( 0 length ( Ht ))) nn = length ( B ) for ( i in 2 : ( nn 1 )) { B [ i ] = 1 / lamda [ i ] + Ht [ i 1 ] + Ht [ i ] } B [ 1 ] = 1 B [ nn ] = 1 B [ fur.loc ] = 1 return ( B ) } BB.I < cmpfun ( BB.Irrig ) #### CC.Irrig < function( Ht ) { fur.loc < FurLoc ( pars constants ) C < as.vector ( rep ( 0 length ( Ht ))) nn = length ( C ) for ( i in 2 : ( nn 1 )) { C [ i ] = Ht [ i ] } C [ 1 ] = 0 C [ nn ] = 0 C [ fur.loc ] = 0 return ( C ) } PAGE 124 124 CC.I < cmpfun ( CC.Irrig ) #### RR.Irrig < function( pars BC1.I BC2.I BC3.I Ht dp fp lamda ET Rf ) { e = ET Re = min ( Rf Rf exp ( pars $ br ( pars $ Tr Ht ))) fur.loc < FurLoc ( pars constants ) R < as.vector ( rep ( 0 length ( Ht ))) for ( i in 1 : ( length ( R ) 1 )) { R [ i ] = ( Ht [ i ] delt e / dp [ i ] + delt Re / fp [ i ] ) / lamda [ i ] } R [ 1 ] = BC1.I R [ length ( R ) ] = BC2.I R [ fur.loc ] = BC3.I return ( R ) } RR.I < cmpfun ( RR.Irrig ) ### Subroutine for Irrigation Phase Iterate.Irrig < function( pars constants BC1.I BC2.I BC3.I OldtHt ET Rf ) { OldmHt = OldtHt if ( Rf > 0 && Rf > ET ) { fp = fPor.Rf ( pars constants OldmHt ( Rf ET )) dp = dPor.Rf ( pars constants OldmHt ( Rf ET )) } else { fp = fPor.ET ( pars constants OldmHt ( ET Rf )) dp = dPor.ET ( pars constants OldmHt ( ET Rf )) } lamda < pars $ Ksat delt / ( fp delz ^ 2 ) A < AA.I ( OldmHt ) B < BB.I ( lamda OldmHt ) C < CC.I ( OldmHt ) R < RR.I ( pars BC1.I BC2.I BC3.I OldtHt dp fp lamda ET Rf ) NewmHt < Tridiag ( A B C R ) return ( NewmHt ) } Iter.Irrig < cmpfun ( Iterate.Irrig ) ##### Subroutine for Drainage Iterate.Drain < function( pars constants BC1.D BC2.D BC3.D OldtHt ET Rf ) { OldmHt = OldtHt fp = fPor.Rf ( pars constants OldmHt ( Rf ET )) dp = dPor.Rf ( pars constants OldmHt ( Rf ET )) if ( Rf == 0 && ET == 0 ) { dp = dPor.ET ( pars constants OldmHt ET ) fp = fPor.ET ( pars constants OldmHt ET ) } if ( Rf > 0 ) { fp = fPor.Rf ( pars constants OldmHt ( Rf )) dp = dPor.Rf ( pars constants OldmHt ( Rf )) } else { fp = fPor.ET ( pars constants OldmHt ( ET )) PAGE 125 125 } lamda < pars $ Ksat delt / ( dp delz ^ 2 ) mHt = OldmHt [ ( L + 1 ) / 2 ] if ( mHt > BC3.D ) { A < AA.D2 ( OldmHt ) B < BB.D2 ( lamda OldmHt ) C < CC.D2 ( OldmHt ) R < RR.D2 ( pars BC1.D BC2.D BC3.D OldtHt dp fp lamda ET Rf ) } else { A < AA.D ( OldmHt ) B < BB.D ( lamda OldmHt ) C < CC.D ( OldmHt ) R < RR.D ( pars BC1.D BC2.D OldtHt dp fp lamda ET Rf ) } NewmHt < Tridiag ( A B C R ) return ( NewmHt ) } Iter.Drain < cmpfun ( Iterate.Drain ) ############## Define Initial Conditions Ht.Ini = RL constants $ Wt.Ini Ht.data < matrix ( 0 N L ) Ht.data [ 1 ] = Ht.Ini BC1.D = RL ( pars $ dwl rowht ) BC2.D = RL ( pars $ dwl rowht ) BC3.D = RL ( pars $ BC3.D rowht ) BC1.I = RL ( pars $ dwl rowht ) BC2.I = RL ( pars $ dwl rowht ) BC3.I = RL ( pars $ BC3.I rowht ) subsur.flux < rep ( 0 N ) pb < tkProgressBar ( title = "progress bar", min = 0 max = N width = 300 ) for ( i in 2 : N ) { if ( Irrig.Data $ Status [ i ] == 0 ) { Ht.data [ i 1 ] = BC1.D Ht.data [ i L ] = BC2.D BC3.D = BC3.D Ht.data [ i ] < Iter.Drain ( pars constants BC1.D BC2.D BC3.D OldtHt = Ht.data [ i 1 ] ET [ i 1 ] Rf [ i 1 ] ) Ht.data [ i ] < ifelse ( Ht.data [ i ] >= ( RL 20 ), ( RL 20 ), Ht.data [ i ] ) subsur.flux [ i ] = flux.DF ( Ht.data [ i L 1 ] ( BC1.D )) } if ( Irrig.Data $ Status [ i ] == 1 ) { Ht.data [ i 1 ] = BC1.I Ht.data [ i L ] = BC2.I BC3.I = BC3.I Ht.data [ i ] < Iter.Irrig ( pars constants BC1.I BC2.I BC3.I Ht.data [ i 1 ] ET [ i 1 ] Rf [ i 1 ] ) subsur.flux [ i ] = flux.DF ( Ht.data [ i L 1 ] ( BC1.I )) PAGE 126 126 } print ( i ) setTkProgressBar ( pb i title = paste ( round ( i / N 100 0 ), "% Progress")) } close ( pb ) Wt.data < RL Ht.data # write.table(Wt.data, file = 'Wtdata.out', sep = ',', col.names = NA, # qmethod = 'double') ## For Writing the simulation result to a output # file # return(Wt.data[N, (L+1)/2]) ##For mid point WT elevation return ( list ( Wt.data = Wt.data slf = subsur.flux )) } PAGE 127 127 LIST OF REFERENCES Abdul, A. S., and R. W. Gillham (1989), Field studies of the effects of the capillary fringe on streamflow generation, Journal of Hydrology 112(1 2), 1 18, doi:10.1016/0022 1694(89)90177 7. Acharya, S., and R. S. Mylavarapu (2011), Selected soil physical properties and implications on water management system in northeast Florida, Soil Sci., 176(1), doi: 10.1097/SS.0b013e31820317c1. Allen, R. G., L. S. Pereira, D. Raes, and M. Smith (1998), Crop Evapotranspiratio n: Guidelines for Computing Crop Requirements, FAO Irrigation and Drainage Paper 56, Food and Agricultural Organization of the U.N., Rome. Bear, J. (1972), Dynamics of Fluids in Porous Media, Elsevier, New York. Bear, J., and A.H. D. Cheng (2008), Modeling groundwater flow and contaminant transport. Theory and Applications of Transport in Porous Media, 23, doi 0.1007/978 1 4020 6682 5. Bierkens, M. F. P. (1998), Modeling water table fluctuations by means of a stoc hastic differential equation, Water Resour. Res., 34(10), 2485 2499, doi:199810.1029/98WR02298. Bouarfa, S., and D. Zimmer (2000), Water table shapes and drain flow rates in shallow drainage systems, Journal of Hydrology, 235(3 4), 264 275, doi:10.1016/S00 22 1694(00)00280 8. Bouwer, H. (1978), Groundwater Hydrology, McGraw Hill Book Co. Bouwer, H., and R. D. Jackson (1974), Determining soil properties, in Drainage for Agriculture, ASA Monograph 17, chap.23, sect. 10, edited by J. van Schilfgaarde 611 672, American Society of Agronomy, Madison, WI. Brooks, R. H., and A. T. Corey (1964), Properties of Porous Media Affecting Fluid Flow, Journal of the Irrigation and Drainage Division, 92(2), 61 90. Brutsaert, W. (1994), The unit response of groundwater outflo w from a hillslope, Water Resour. Res., 30(10), 2759 2763, doi:10.1029/94WR01396. Brutsaert, W. (2005), Hydrology: An Introduction, 378 379 pp., Cambridge Univ. Press, Cambridge, U.K. Butler, J. J., Jr., G. J. Kluitenberg, D. O. Whittemore, S. P. Loheide I I, W. Jin, M. A. Billinger, and X. Zhan (2007), A field investigation of phreatophyte induced fluctuations in the water table, Water Resour. Res., 43, W02404, doi:10.1029/2005WR004627. Campbell, K. L., J. S. Rogers, and D. R. Hensel (1978), Water table con trol for potatoes in Florida, Transactions of the American Society of Agricultural Engineering, 21, 701 705. PAGE 128 128 Carsel, R., and R. Parrish (1988), Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24(5), 755 769, doi:10.1029/WR024i005p00755. Cartwright, N., P. Nielsen, and P. Perrochet (2005), Influence of capillarity on a simple harmonic oscillating water table: Sand column experiments and modeling, Water Resour. Res., 41, W08 416, doi:10.1029/2005WR0040 23. Water Resour. Res., 31 (9), 2377 2378, doi:199510.1029/95WR01488. Chapman, T. G., and R. F. Dressler (1984), Unsteady shallow groundwater f low over a curved impermeable boundary, Water Resour. Res., 20(10), 1427 1434, doi:10.1029/WR020i010p01427. Childs, E. C. (1960), The Nonsteady State of the Water Table in Drained Land, J. Geophys. Res. 65 (2), PP. 780 782, doi:196010.1029/JZ065i002p00780. Childs, E. C. (n.d.), Drainage of Groundwater Resting on a Sloping Bed, Water Resour. Res., 7(5), P. 1256, doi:197110.1029/WR007i005p01256. Chow, V. T. (1964), Handbook of Applied Hydrology McGraw Hill, New York. Dolan, T. J., A. J. Hermann, S. E. Bayley and J. Zoltek Jr. (1984), Evapotranspiration of a Florida, U.S.A., freshwater wetland, Journal of Hydrology 74 (3 4), 355 371, doi:10.1016/0022 1694(84)90024 6. Duke, H. R. (1972), Capillary properties of soils influence upon specific yield, Trans. ASAE 688 691, FAWN, (2011), Florida Automated Weather Network, UF/IFAS Extension, University of Florida, Gainesville, Florida, USA. Feddes, R. A., E. Bresler, and S. P. Neuman (1974), Field test of a modified numerical model for water uptake by root systems, W ater Resour. Res., 10(6),1199 1206, doi:197410.1029/WR010i006p01199. Feddes, R. A., E. Bresler, and S. P. Neuman (1976), Field test of a modified numerical model for water uptake by root systems, Water Resour. Res., 10 (6), 1199 1206, doi:197410.1029/WR010i 006p01199. Feddes, R.A., P.J. Kowalik, and H. Zaradny (1978), Simulation of Field Water Use and Crop Yield John Wiley & Sons. New York, Freeze, R. A. (1971), Three Dimensional, Transient, Saturated Unsaturated Flow in a Groundwater Basin, Water Resour. Re s., 7(2), 347 366, doi:10.1029/WR007i002p00347. PAGE 129 129 Freeze, R. A., and Cherry, John A (1979), Groundwater Prentice Hall, Old Tappan, N. J. Gardner, W. R. (1958), Some steady state solutions of the unsaturated moisture flow equation with application to evapor ation from a water table, Soil Sci. 85 (4). Ghezzehei, T. A., T. J. Kneafsey, and G. W. Su (2007), Correspondence of the Gardner and van Genuchten Mualem relative permeability function parameters, Water Resour. Res., 43, W10417, doi:10.1029/2006WR005339. Gillham, R.W. (1984), The capillary fringe and its effect on water table response, Journal of Hydrology 67 (1 4), 307 324, doi:10.1016/0022 1694(84)90248 8. Glover, R. E. (1966), Ground water movement, United States Bureau of Reclamation. Engineering monog raphs, no. 31, U.S. Dept. of the Interior, Bureau of Reclamation. Gribovszki, Z., J. Szilgyi, and P. Kalicz (2010), Diurnal fluctuations in shallow groundwater levels and streamflow rates and their interpretation A review, Journal of Hydrology 385 (1 4), 371 383, doi:10.1016/j.jhydrol.2010.02.001. Gribovszki, Z., P. Kalicz, J. Szilgyi, and M. Kucsara (2007), Riparian zone evapotranspiration estimation from diurnal groundwater level fluctuations, Journal of Hydrology 349 (1 2), 6 17, doi:10.1016/j.j hydrol.2007.10.049. Harbaugh, A.W. (2005), MODFLOW 2005, The U.S. Geological Survey modular ground water model the ground water flow Process: U.S. Geological Survey Techniques and Methods 6 A16, 6.16 6.18. Healy, R., and P. Cook (2002), Using groundwater levels to estimate recharge, Hydrogeol J 10 (1), 91 109, doi:10.1007/s10040 001 0178 0. by Nick Cartwright et al., Water Resour. Res., 42, W11601, doi:10.1029/2006WR005042. Hilberts, A. G. J., E. E. van Loon, P. A. Troch, and C. Paniconi (2004), The hillslope storage Boussinesq model for non constant bedrock slope, Journal of Hydrology 291 ( 3 4), 160 173, doi:10.1016/j.jhydrol.2003.12.043. Hilberts, A. G. J., P. A. Troch, and C. Paniconi (2005), Storage dependent drainable porosity for complex hillslopes, Water Resour. Res., 41 13 PP., doi:200510.1029/2004WR003725. Hilberts, A. G. J., P. A. Troch, C. Paniconi, and J. Boll (2007), Low dimensional modeling of hillslope subsurface flow: Relationship between rainfall, recharge, and unsaturated storage dynamics, Water Resour. Res., 43, W03445, doi:10.1029/2006WR004964. Hillel, D. (1998), Environme ntal Soil Physics Academic Press, San Diego. PAGE 130 130 Jones, J.W., L.H. Allen, S.F. Shih, J.S. Rogers, L.C. Mammond, A.G. Smajstrala, and J.D. Martsolf. 1984. Estimated and measured evapotranspiration for Florida Climate, Crops and Soils. Agric. Exp. Stn., IFAS, U niv. of Florida, Gainesville. Khan, S., and K. R. Rushton (1996a), Reappraisal of flow to tile drains II. Time variant response, Journal of Hydrology, 183(3 4), 367 382, doi:10.1016/0022 1694(95)02948 6. Khan, S., and K. R. Rushton (1996b), Reappraisal of flow to tile drains III. Drains with limited flow capacity, Journal of Hydrology, 183(3 4), 383 395, doi:10.1016/0022 1694(95)02949 4. Lai, C.T., and G. Katul (2000), The dynamic role of root water uptake in coupling potential to actual transpiration, Adv. Water Resour 23 (4), 427 439, doi:10.1016/S0309 1708(99)00023 8. Iturbe (2009), Ecohydrology of groundwater dependent ecosystems: 1. Stochastic water table dynamics, Water Resour. Res., 45, W05419, doi:10.1029/2008WR007292. Lautz, L. K. (2007), Estimating groundwater evapotranspiration rates using diurnal water table fluctuations in a semi arid riparian zone, Hydrogeol J 16 (3), 483 497, doi:10.1007/s10040 007 0239 0. Loheide II, S. P. (2008 ), A method for estimating subdaily evapotranspiration of shallow groundwater using diurnal water table fluctuations, Ecohydrol. 1 (1), 59 66. Loheide II, S. P., J. J. B. Jr, and S. M. Gorelick (2005), Estimation of groundwater consumption by phreatophytes using diurnal water table fluctuations: A saturated unsaturated flow assessment, Water Resour. Res., 41, 14 PP., doi:200510.1029/2005WR003942. Lowry, C. S., and S. P. L. II (2010), Groundwater dependent vegetation: Quantifying the groundwater subsidy, Wat er Resour. Res., 46 8 PP., doi:201010.1029/2009WR008874. Luo, Y., and M. Sophocleous (2010), Seasonal groundwater contribution to crop water use assessed with lysimeter observations and model simulations, Journal of Hydrology 389 (3 4), 325 335, doi:10.10 16/j.jhydrol.2010.06.011. Milly, P. C. D. (1987), Estimation of Brooks Corey Parameters from water retention data, Water Resour. Res., 23(6), 1085 1089, doi:10.1029/WR023i006p01085. Monteith, J. L. (1965), Evaporation and environment, Symp. Soc. Exp. Biol. 19 205 234. Moody, W. T. (1966), Nonlinear Differential Equation of Drain Spacing, Journal of the Irrigation and Drainage Division 92 (2), 1 10. PAGE 131 131 Morel Seytoux, H. J., P. D. Meyer, M. Nachabe, J. Tourna, M. T. van Genuchten, and R. J. Lenhard (1996), Par ameter Equivalence for the Brooks Corey and Van Genuchten Soil Characteristics: Preserving the Effective Capillary Drive, Water Resour. Res., 32(5), 1251 1258, doi:10.1029/96WR00069. Mould, D. J., E. Frahm, T. Salzmann, K. Miegel, and M. C. Acreman (2010) Evaluating the use of diurnal groundwater fluctuations for estimating evapotranspiration in wetland environments: case studies in southeast England and northeast Germany, Ecohydrol 3 (3), 294 305, doi:10.1002/eco.108. Mualem, Y. (1984), A modified depen dent domain theory of hysteresis, Soil Sci., 137(5), 283 291. Naumburg, E., R. Mata gonzalez, R. G. Hunter, T. Mclendon, and D. W. Martin (2005), Phreatophytic Vegetation and Groundwater Fluctuations: A Review of Current Research and Application of Ecosyst em Response Modeling with an Emphasis on Great Basin Vegetation, Environmental Management, 35(6), 726 740, doi:10.1007/s00267 004 0194 7. Munoz Arboleda, F., R. Mylavarapu, C. Hutchinson, and K. Portier (2008), Nitrate nitrogen concentrations in the perche d ground water under seepage irrigated potato cropping systems, J. Environ.Qual. 37 (2), 387 394, doi:10.2134/jeq2006.0545. Nachabe, M. H. (2002), Analytical expressions for transient specific yield and shallow water table drainage, Water Resour. Res., 38( 10), 1193, doi:10.1029/2001WR001071. Nachabe, M., C. Masek, and J. Obeysekera (2004), Observations and Modeling of Profile Soil Water Storage above a Shallow Water Table, Soil Sci. Soc. Am. J. 68 (3), 719 724. Nachabe, M., N. Shah, M. Ross, and J. Vomacka (2005), Evapotranspiration of two vegetation covers in a shallow water table environment, Soil Sci. Soc. Am. J. 69, 492, doi:10.2136/sssaj2005.0492. Nash, J. E., and J. V. Sutcliffe (1972), River flow forecast ing through conceptual models part I: A discussion of principles, Journal of Hydrology, 10(3), 282 290, doi:10.1016/0022 1694(70)90255 6. Neuman, S. P. (1972), Theory of flow in unconfined aquifers considering delayed response of the water table, Water Res our. Res., 8 (4), doi:197210.1029/WR008i004p01031. Neuman, S. P. (1987), On methods of determining specific yield, Ground Water 25 (6), 679 684, doi:10.1111/j.1745 6584.1987.tb02208.x. Nichols, W. D. (1993), Estimating discharge of shallow groundwater by tr anspiration from greasewood in the Northern Great Basin, Water Resour. Res., 29(8), 2771 2778, doi:10.1029/93WR00930. PAGE 132 132 Nichols, W. D. (1994), Groundwater discharge by phreatophyte shrubs in the Great Basin as related to depth to groundwater, Water Resour. R es., 30(12), 3265 3274, doi:10.1029/94WR02274 Nielsen, P., and P. Perrochet (2000), Water table dynamics under capillary fringes: Experiments and modelling, Adv. Water Resour., 23(1), 503 515. doi:10.1016/S0309 1708(99)00038 X Novakowski, K. S., and R. W. Gillham (1988), Field investigations of the nature of water table response to precipitation in shallow water table environments, Journal of Hydrology 97 (1 2), 23 32, doi:10.1016/0022 1694(88)90063 7. NRCS, (1999), Soil Taxonomy United States Department of Agriculture US Government Printing Office, Washington DC. Pandey, R. S., A. K. Bhattacharya, O. P. Singh, and S. K. Gupta (1997), Water table drawdown during drainag e with evaporation/evapotranspiration, Agric Water Manage 35(1 2), 61 73, doi:10.1016/S0378 3774(97)00032 2. Park, E., and J. C. Parker (2008), A simple model for water table fluctuations in response to precipitation, Journal of Hydrology 356 (3 4), 344 349, doi:10.1016/j.jhydrol.2008.04.022. Philip, J. R. (1969), Theory of infiltration, Adv. Hydrosci., 5, 215 296 Philip, J. R. (1985), Approximate analysis of the borehole permeameter in unsaturated soil, Water Resour. Res., 21 (7), 1025 1033, doi:198510.1 029/WR021i007p01025. Pikul, M. F., R. L. Street, and I. Remson (1974), A numerical model based on coupled one dimensional Richards and Boussinesq equations, Water Resour. Res., 10(2), 295 302, doi:10.1029/WR010i002p00295. Pitts, D.J., Smajstrla A.G., 1989. Irrigation systems for crop production in Florida: descriptions and costs. Circular 821. Fl. Coop. Ext. Serv. Gainesville, Fl. Priestley, C. H. B., and R. J. Taylor (1972), On the assessment of surface heat flux and evaporation using large scale parameter s, Mon. Weather Rev. 100 81 92 Pumo, D., S. Tamea, L. V. Noto, F. Miralles Wilhem, and I. Rodriguez Iturbe (2010), Modeling belowground water table fluctuations in the Everglades, Water Resour. Res., 46, W11557, doi:10.1029/2009WR008911. R Development Core Team (2011), R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3 900051 07 0, URL http://www.R project.org. Raats, P. A. C. ( 1970), Steady infiltration from line sources and furrows, Soil Sci. Soc. Am. Proc., 34, 709 714. PAGE 133 133 Reynolds, W. D., and D. E. Elrick (1985), In situ measurements of field saturated hydraulic parameter using Guelph Permeam eter, Soil Sci., 140 (4), 292 302. Reynolds, W. D., and D. E. Elrick (1987), A laboratory and numerical assessment of the guelph permeameter method, Soil Sci., 144(4), 282 299. Ritchie, J. T. (1972), Model for predicting evaporation from a row crop with i ncomplete cover, Water Resour. Res., 8(5), 1204 1213, doi:10.1029/WR008i005p01204. Ritzema, H. P. (1994), Drainage principles and applications International Institute for Land Reclamation and Improvement, Netherlands. Rosa, J.A. (2000), Modeling water ta ble response to Subirrigation with a buried microirrigation line source for potato production, PhD Dissertation, Dep. of Ag. and Bio. Eng. Univ. of Florida, Gainesville, Florida, USA. Rucker, D. F., A. W. Warrick, and T. P. A. Ferr (2005), Parameter equiv alence for the Gardner and van Genuchten soil hydraulic conductivity functions for steady vertical flow with inclusions, Adv. in Water Resour. 28 (7), 689 699, doi:10.1016/j.advwatres.2005.01.004. Rupp, D. E. and J. S. Selker (2006), On the use of the Bous sinesq equation for interpreting recession hydrographs from sloping aquifers, Water Resour. Res., 42, W12421, doi:10.1029/2006WR005080. Russo, D. (1988), Determining soil hydraulic properties by parameter estimation: On the selection of a model for the hyd raulic properties, Water Resour. Res., 24(3), 453 459, doi:10.1029/WR024i003p00453. Schilling, K. E., and J. R. Kiniry (2007), Estimation of evapotranspiration by reed canarygrass using field observations and model simulations, Journal of Hydrology 337 (3 4), 356 363, doi:10.1016/j.jhydrol.2007.02.003. Shah, N., M. Nachabe, and M. Ross (2007), Extinction depth and evapotranspiration from ground water under selected land covers. Ground Water, 45, 329 338. doi:10.1111/j.1745 6584.2007.00302.x. Simunek, J., M. Sejna, and M.Th. van Genuchten (1998), The HYDRUS 1D: Software package for simulating the one dimensional movement of water, heat, and multiple solutes in variably saturated media.Version 2.0. Riverside, California: US Salinity Laboratory, ARS/USDA. Singh, S., and C. Jaiswal (2006), Numerical Solution of 2D Free Surface to Ditch Drains in Presence of Transient Recharge and Depth Dependent ET in Sloping Aquifer, Water Resources Management, 20(5), 779 793, doi:10.1007/s11269 005 9007 x. PAGE 134 134 Singl eton, V. D. 1990. Investigation of potato water use in the Tri county agricultural area of Putnam, St. Johns, and Flagler Counties, Florida. SJ 90 13.Tech. Publ. Johns River Water Management District. Palatka, FL. Skaggs, R. (1973), Water table movement du ring subirrigation, Trans. of ASAE, 16, 988 993. Skaggs, R. W. (1980), A water management model for shallow water table soils. Technical Report 134 University of North Carolina, Raleigh, NC. Skaggs, R. W. (1991), Modeling water table response to subirrig ation and drainage, Trans. of ASAE, 34(1) 169 175. Smajstrla, A., S. Locascio, D. Weingartner, and D. Hensel (2000), Subsurface drip irrigation for water table control and potato production, Appl. Eng. Agric. 16 (3), 225 229. Sophocleous, M. A. (1987), Th e Role of Specific Yield in Ground Water Recharge Estimations: A Numerical Study, Ground Water 23 (1), 52 58. Sophocleous, M. A. (1991), Combining the soil water balance and water level fluctuation methods to estimate natural groundwater recharge: Practica l aspects, Journal of Hydrology, 124(3 4), 229 241, doi:10.1016/0022 1694(91)90016 B. Soylu, M. E., E. Istanbulluoglu, J. D. Lenters, and T. Wang (2010), Quantifying the impact of groundwater depth on evapotranspiration in a semi arid grassland region, Hydrol. Earth Syst. Sci. 7 6887 6923, doi:10.5194/hessd 7 6887 2010. Stauffer, F., H. J. Franke, and T. Dracos (1992), Hysteretic storativity concept for aquifer simulation, Water Resour. Res., 28 (9),2307 2314, doi:199210.1029/92WR01148. Tamea, S., F. La Iturbe (2009), Ecohydrology of groundwater dependent ecosystems: 2. Stochastic soil moisture dynamics, Water Resour. Res., 45, W05420, doi:10.1029/2008WR007293. Tang, Y. K. and R. W. Skaggs (1977), Experiment al evaluation of theoretical solutions for subsurface drainage and irrigation, Water Resour. Res., 13(6), 957 965, doi:10.1029/WR013i006p00957. Tritscher, P., W. W. Read, and P. Broadbridge (2000), Specific yield for a two dimensional flow, Water Resour. Re s., 36 (6), PP. 1393 1402, doi:200010.1029/1999WR900333. Troch, P. (1992), Conceptual basin scale runoff process models for humid catchments: Analysis, synthesis and applications Ph.D. Thesis p 62 63. Ghent Univ., Ghent, Belgium. Troch, P. A., C. Paniconi and E. E. van Loon (2003), Hillslope storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 1. Formulation and characteristic response, Water Resour. Res., 39, 12 PP., doi:200310.1029/2002WR001728. PAGE 135 135 Troch, P. A., F. P. De Troch, and W. Brutsaert (1993), Effective water table depth to describe initial conditions prior to storm rainfall in humid regions, Water Resour. Res., 29(2), 427 434, doi:10.1029/92WR02087. Van Beers, W.F.J. 1979. Some nomographs for the calcula tion of drain spacing 3rd ed. ILRI Bulletin 8, Wageningen, 46 p. van Dam, J. C., and R. A. Feddes (2000), Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation, Journal of Hydrology, 233(1 4), 72 85, d oi:10.1016/S0022 1694(00)00227 4. van Dam, J. C., J. H. M. Wsten, and A. Nemes (1996), Unsaturated soil water movement in hysteretic and water repellent field soils, Journal of Hydrology, 184(3 4), 153 173, doi:10.1016/0022 1694(95)02996 6. van der Molen, W. H., and J. Wesseling (1991), A solution in closed form and a series solution formula, Agricultural Water Management 19 (1), 1 16, doi:10.1016/0378 3774(91)9005 8 Q. van Genuchten, M. T. (1980), A closed form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (5), 892, doi:10.2136/sssaj1980.03615995004400050002x. van Schilfgaarde, J. (1974) Nonsteady flow to drains, I n Drainage for Agriculture, ASA Monograph 17, chap.23, sect. 10, edited by J. van Schilfgaarde, 611 672, American Society of Agronomy, Madison, WI. Verhoest, N. E. C., and P. A. Troch (2000), Some analytical solutions of the linearized Boussinesq equation with recharge for a sloping aquifer, Water Resour. Res., 36(3), P. 793, doi:200010.1029/1999WR900317. Verhoest, N. E. C., V. R. N. Pauwels, P. A. Troch, and F. P. De Troch (2002), Analytical Solution for Transient Water Table Heights and Outflows from Inc lined Ditch Drained Terrains, Journal of Irrigation and Drainage Engineering, 128(6), 358 364, doi:10.1061/(ASCE)0733 9437(2002)128:6(358). Warrick, A. W. (1995), Correspondence of hydraulic functions for unsaturated soils, Soil Sci. Soc. of Am. J. 59 (2), 292, doi:10.2136/sssaj1995.03615995005900020003x. Weir, G. J. (1987), Steady infiltration from small shallow circular ponds, Water Resour. Res., 23(4), 733 736, doi:10.1029/WR023i004p00733. White, W.N. 1932. A method of estimating ground water supplies b ased on discharge by plants and evaporation for soil. USGS Water Supply Paper 659: 1 105. PAGE 136 136 Woods, R. A., M. Sivapalan, and J. S. Robinson (1997), Modeling the spatial variability of subsurface runoff using a topographic index, Water Resour. Res., 33(5), 106 1 1073, doi:10.1029/97WR00232 Zhu, J., and B. P. Mohanty (2004), Soil hydraulic parameter upscaling for steady state flow with root water uptake, Vadose Zone Journal 3 (4), 1464 1470, doi:10.2113/3.4.1464. Zhu, J., and B. P. Mohanty (2002), Upscaling of soil hydraulic properties for steady state evaporation and infiltration, Water Resour. Res., 38(9), 1178, doi:10.1029/2001WR000704. PAGE 137 137 BIOGRAPHICAL SKETCH Subodh Acharya was born in 1980 in Nepal, in a small village called Shankhu, approximately 36 kilometers from the capital city Kathmandu. He earned his b degree in Agriculture from the Tribhuvan University, Nepal In fall 2006, he joined the graduate program at the University of Florida, Soil and Water Science De partment (SWSD) as a m s tudent. After earning his m in soil and water science in summer 2008, h e joined the PhD prog ram in the SWSD department in fall 2008 and received the Graudate Alumni Fellowship to support his study This thesis i s the culmination of the studies h e conducted during 2008 2011 as part of his PhD dissertation research His primary research interests include the study of soil water plant interactions in areas with shallow water table through field as well as modeling b ased approaches. He is also highly interested in the surface and subsurface hydrology, especially in shallow water table areas, both natural and managed. In May 20 12 h e successfully defended his dissertation and earned the doctorate degree in soil and wat er science in August 2012 